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SUMMARY 


A numerical study of non-reactive and reactive 
axisymmetric combustor flows with and without swirl is 
presented. Closure of the Reynolds equations is achieved by 
three different levels of models: k-e, algebraic stress and 
Reynolds stress closure. Performance of two locally non- 
equilibrium and one equilibrium algebraic stress models is 
analyzed assuming four different pressure-strain models. A 
comparison is also made of the performance of a high and a 
low Reynolds number model for combustor flow calculations 
using Reynolds stress closures. Effects of diffusion and 
pressure-strain models on these closures are also 
investigated. Two different models for the scalar transport 
are presented. One employes the second-moment closure which 
solves the transport equations for the scalar fluxes, while 
the other solves the algebraic equations for the scalar 
fluxes. In addition, two cases of non-premixed and one case 
of premixed combustion are considered. Fast- and finite-rate 
chemistry models are applied to non-premixed combustion. 
Both models show promise for application in gas turbine 
combustors. However, finite rate chemistry models, which are 
more realistic, need to be further examined to establish a 
suitable coupling of the heat release effects on the 
turbulence field and the rate constants. 


xii 


CHAPTER 1 


INTRODUCTION 


1.1 BACKGROUND 

The calculation of turbulent combusting flows has 
received considerable attention in recent years. This is due 
to different reasons. Some of these are: demands for higher 
efficiency combustors, higher cost of measurement 
instruments, rapid growth of computer technology, limitation 
of fossil fuel resources, and a better understanding of 
pollutant formations which jeopardize the ecology in 
industrial societies. As a result, an increased interest in 
combustion modelling has been generated by researchers in 
this area. 

At first, attention is only given to the handling of 
the complexities of real chemical reactions, which involve 
scores of individual species (Lavoie, et al.,1970). In 
reality, fifty to one hundred elementary reactions, 
involving ten to fifty chemical species, are required for a 
complete description of the kinetic process. Typical 
examples are concerned with the accurate calculation of the 
equilibrium composition of the gases in a rocket nozzle, 
followed by the kinetically-influenced variation of the 
composition of the gases as they flow through the nozzle. 


Later, two- and three-dimensional variations in space and 
variations in time are taken into account (Butler and 
O'Rouke, 1977; Griffin et al. f 1978), and the development of 
"turbulence models" permits realistic representation of the 
flow patterns and chemical reactions. 

In order to construct a comprehensive combustion model, 
the problem may be broadly grouped into three categories: 
numerical methods, turbulence models, and chemical kinetics. 
There are some serious numerical obstacles involved in a 
comprehensive combustion model. Time scales associated with 
molecular transport and chemical reactions are vastly 
different, and there are sharp gradients in the temperature 
and density fields. Simulation of such a process would 
involve a set of partial differential equations that include 
mass, momentum, energy and species conservation. This set of 
equations describing the combustion flow field are coupled 
and generally they cannot be easily solved by conventional 
procedures because of extremely disparate physical scales 
introduced by steep gradients involved in the reaction 
processes. For example, in the work done by Dwyer and 
Sanders (1977), a numerical technique known as operator 
splitting was used to compensate for the vastly different 
time scales associated with molecular transport, chemical 
transport and unsteady wave propagation. The fluid dynamics 
portion of the model problem was solved using explicit 
finite difference scheme, while the split reaction terms 
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were solved by an ordinary differential equation technique. 
Otey (1978) has presented a review of different numerical 
methods. His work retained the essential features of 
chemical flows and allowed insight into the numerical 
difficulties involved in the solution of such a complex 
system. 

However, among the most important questions in current 
combustion research are those regarding the role of 
interaction between the fluid mechanics and the chemical 
reactions (Smoot and Hill, 1983). Chemical reactions take 
place in the molecular level. On the other hand, turbulent 
motion plays an important role in the mixing of the reacting 
species. Local turbulence controls the time and frequency 
that each of the reactants and products are mixed together, 
thus allowing reaction to proceed. In turn, the local 
instantaneous reaction processes themselves, often 
associated with local heat release or absorption, density 
and volume changes, etc. can impact the local turbulent fluid 
mechanics . 

In practice, the fluid mechanics and chemical kinetics 
are so complex that various modelling assumptions are 
required to render the solution of the conservation 
equations possible. These assumptions which combine physical 
reality with ease of formulation are concerned with the 
turbulent nature of the flow, the flame, the combustion 
characteristics, and the radiative heat transfer from the 
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products of combustions. Improper models result in erroneous 
prediction of combustion efficiency, temperature 
distribution in the combustion system and pollutant 
formations . 

Numerous turbulence models have successfully been 
devised for constant properties flows (Launder et al.,1975; 
Lumley, 1975a; Reynolds , 1976) . These provide the basis for 
extension to flows with variable properties and combustion. 
However, one of the main problem in devising physically 
valid models for reactive flows lies in the characterization 
of the time-mean net rate of formation or destruction of the 
molecular species due to chemical reactions (Pratt, 1979 
;Bilger, 1980) . Although the kinetic mechanisms are not 
always known and kinetic rate constants are difficult to 
identify, the major problem lies not in these areas but in 
obtaining the proper time-mean rate of formation or 
destruction of molecular species due to the presence of 
turbulence (Borghi, 1974) . Analytical expressions relating 
the instantaneous reaction rates to associated flow 
quantities are always highly non-linear functions of the 
temperature and species concentrations, and thus the 
knowledge of the mean values of these latter quantities is 
insufficient to allow the evaluation of the mean formation 
rates (Pratt , 1979) . In order to express the mean formation 
rate in a correct manner, it is convenient to use the 
technique of Reynolds decomposition and express the 
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Arrhenius term in terms of an infinite series involving the 
mean values of the scalars and their higher statistics and 
cross-correlations. This series is not always convergent 
because the correlation terms can be of comparable magnitude 
with the mean quantities (Borghi, 1974) . Depending on the 
type of reactions considered and on the chemical species 
involved, Borghi (1974) estimated that at least seven terms, 
involving third and lower order moments of the scalars, in 
the expansion have to be considered. For any realistic 
combustion scheme, many equations of this type are present. 
Therefore, it is clear that closure of turbulent reacting 
flow is a formidable problem. This problem was recognized by 
various researchers and different alternate paths have been 
suggested. A review of the state of the art of different 
combustion models for turbulent diffusion flames has been 
given by Bilger (1976) . 

The effects of these chemistry-turbulence interactions 
are quite different for different types of chemical 
reactions. To help identify these effects, it is convenient 
to identify two hypothetical time scales: the reaction time 
scale and the turbulent time scale. The reaction time scale 
is defined as a typical time for the reacting species of 
interest to react completely to its equilibrium value. The 
turbulence time scale is chosen to be a typical fine-scale 
mixing time for scale reduction by turbulent breakup of 
large eddies. This time scale must be adequate for molecular 
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interaction to take place (micromixing) . The turbulence time 
scale then is the time required for mixing to proceed to the 
molecular level before reaction can occur. The reaction time 
scale is the time required for these reacting species, once 
contacted, to react completely to form their products. 
Approaches for incorporating chemical reactions in turbulent 
systems can be characterized by examining the relationship 
between these two time scales. 

If the reaction time scale, t p , is much greater than 
the turbulent time scale, t t , then the reactions are very 
slow compared to changes in the local turbulence. In this 
case, if the fluctuations of any variable are very small, 
the chemistry is unaware of the presence of the turbulent 
fluid mechanics, and the effect of the fluctuations on the 
reaction rate can be ignored. However, the mean reaction 
rate is highly sensitive to temperature fluctuations. Thus, 
even though reaction rates are slow but significant 
temperature fluctuations exist, still the mean reaction rate 
should account for turbulence effects. Therefore, only in 
the very special case is the mean reaction rate equal to the 
reaction rate calculated from the mean variables. Although 
this approximation has been used by many researchers, it has 
been shown to be valid only in a very limited number of 
cases. When reaction rates are not suffeciently slow 
compared to the local turbulence, the use of this 
approximation produces appreciable error (Pratt, 1979) . 
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When the reaction time scale is of the same order of 
magnitude as the turbulent time scale, both the chemical 
kinetics and the turbulent fluctuations must be accounted 
for. It is this area that has been identified by combustion 
researchers as the area that has the most need of specific 
research advances (Smoot and Hill, 1983). Recently, Bilger 
(1978) reported an important experimental observation 
relating to chemistry/turbulence interaction in a diffusion 
flames. Over a broad region of a given non-equilibrium, 
laminar, hydrocarbon diffusion flame, the molecular species 
composition was only a function of the equivalence ratio, 
eVen though the products were not in thermodynamic 
equilibrium. He concluded that the species reaction rate 
seems to be only a function of the mixture fraction in these 
flames. This observation holds great significance for the 
interactions between kinetically limited chemistry and the 
turbulent flow field. If the expeimental data base could be 
extended to turbulent flames with similar relationships 
between local instantaneous chemical composition and the 
local instantaneous conserved scalar, then statistical 
techniques for calculating mean composition could be 
extended to non-equilibrium combustion (Bilger, 1980) . 

However, reactions associated with heat release in the 
high temperature oxidation of hydrocarbon fuels have time 
scales very short compared to the time scale of the 
turbulent micromixing process (Libby and Williams, 1980) . In 
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this case, the reactions occur quickly once the reacting 
species are mixed together. The fast chemistry assumption 
can then be applied. Unfortunately this assumption can be 
applied only to situatioins where the fuel and oxidant enter 
in separate streams. In this type of turbulent reacting 
flows, the assumption can be made that the micromixing 
process is rate limiting and not the chemical kinetic 
process. Thus, as far as the overall flame is concerned, the 
chemistry is fast enough to be considered in local 
instantaneous equilibrium. Therefore, fuel and oxidant 
cannot both exist at the same point. Once the reactants are 
mixed through turbulent motions, the reactions proceed to 
equilibrium instantaneously. If the assumptions are made 
that all species and heat diffuse at the same rate and the 
heat loss to the surrounding is negligible compared with the 
heat release, then the instantaneous chemical composition 
and temperature can be determined in terms of a single 
conserved scalar quantity. For these cases, the conventional 
conserved scalar or mixture fraction is defined to identify 
the degree of "mixedness" at a point. With these 
simplifications, the reacting flow problem is reduced to an 
equivalent non-reacting mixing problem (Bilger, 1980) . A 
significant weakness in this approach is that no details 
concerning the formation and emission of carbon monoxide, 
nitric oxide, and unburnt fuel are available. All these 
require consideration of finite rate chemistry. However, 
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this method is sufficient even when equilibrium cannot be 
assumed, but the reaction rates are sufficiently fast that 
chemical composition is only a function of the local 
equivalance ratio, and thus only a function of the mixture 
fraction. Physically, this implies that the reactions 
proceed quickly to some condition that approximates the 
equilibrium condition (Smith and Smoot, 1981) . 

The mean density can be evaluated once the scalar 
probability density function (PDF) is known. With 

appropriate incorporation of the intermittency of pure fuel 
or oxidant streams, some effects of turbulence on the 
chemical reactions can be accounted for. The approach 
adopted involves specifying a two-parameter PDF. The mean 
and variance of the conserved scalar which are determined 
from solutions of their respective transport equations can 
be used to determine the unknowm parameters. Different two- 
parameter PDF's have been proposed and tested by a number of 
researchers (Spalding, 1971 ; Lockwood and Naguib, 1975 ; 

Rhode, 1975). However, evidence in support of the beta- 

distribution has been provided by Jones (1977) in his 
calculation of diffusion flames. 

Despite the rapid advances made in computer technology, 
a direct solution of the time-dependent conservation 
equations for turbulent flows is not currently practical. 
For example, the calculation of the diurnal cycle by 
Deardorf (1974) using a sub-grid-scale scheme required a 
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week's computing time using the whole resources of a CDC 
7600 computer. The conventional technique is to solve the 
time-averaged equations. In this way, the equations become 
identical to the instantaneous form of the equations only in 
the time-mean variable, but there are a large number of 
extra terms involving the fluctuating components (Borghi, 
1974) . Terms such as pu^Uj ,Up u^, etc. appear. A constant 
density flow simplifies these problems greatly since p'=0. 
For variable density flow, more often, it is assumed that p' 
is not correlated with u^, thus, terms involving pare 
neglected and the equations reduce to those for the constant 
density flow. However, some measurements cited by Bilger 
(1976) indicate that terms such as Up v can be of the same 
order as p'u^Uj and sometimes greater than the turbulent 
momentum flux puv . However, there is a way to circumvent 
this difficulty. Bilger (1975; 1976) points out how Favre 
averaging can be applied to deal with this problem. In Favre 
averaging, quantities are weighted by the instantaneous 
density before averaging. The resulting partial differential 
equations are identical in form to the uniform density flow 
equations, except Favre-averaged variables replace the 
conventional Reynolds-averaged ones (Favre, 1969) . 
Furthermore, the density remaining in the equation is still 
the time mean density. 

In the fluid mechanics area, a survey of the mean- 
turbulence-field closure models by Mellor and Herring (1973) 
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gives an excellent discussion of what has been achieved in 
the mathematical modelling of turbulence up to 1973. 
Spalding (1975) gives a discussion of solved and unsolved 
problems in turbulence modelling. He focuses on the k-e 
model and enumerates its advantages and shortcomings. It is 
the deficiencies in the k-€ model which have encouraged 
developments in more advanced turbulence modelling. 
Different classes of turbulence models have been developed 
such as Reynolds stress closures (Hanjalic and 
Launder, 1972) , algebraic stress/scalar flux closures (Mellor 
and Yamada,1974 ; Rodi, 1976; Gibson and Launder, 1976), 
sub-grid-scale scheme (Deardorff, 1973,1975; Schumann, 1975 ; 
Kwak et al., 1975), etc. All of these focus on the non- 
isotropic nature of the eddy viscosity. 

In Reynolds stress closures, the stress components are 
no longer related to the local mean strain rates but are 
determined from their transport equations (Launder, 1979) . A 
similar approach is used to determine the turbulent scalar 
fluxes. Models employing transport equations for the 
individual turbulent stress and flux components simulate the 
turbulent processes most realistically and are therefore 
potentially better than simpler models. However, they are 
not thoroughly tested and computationally more expensive, 
and hence at the present state of development not very 
suitable for practical applications. However, they are 
important as a starting point for deriving algebraic 
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expressions for the turbulent stresses and fluxes 
(Rodi,1976; Mellor and Yamada, 1974) . It seems that such 
expressions, together with k and e equations, are sufficient 
for most engineering problems and that there are hardly any 
cases where a full transport equation model is needed. They 
combine, at least to some extent, the numerical simplicity 
of the k-e model with the generality of the Reynolds stress 
closures when it comes to accounting for the effect of body 
forces such as buoyancy and rotation. Finally, a completely 
different turbulence modelling approach is to solve the 
Navier-Stokes equations directly by simulation and modelling 
only the sub-grid-scale turbulence. The present computers 
are too small and too slow to resolve the small-scale 
turbulent motion in a numerical solution. However, the 
computer capacity is sufficient to solve the time-dependent 
equations for the large scale motions; the small-scale 
turbulence that cannot be resolved with the chosen numerical 
grid must then be approximated by a model. The small-scale 
turbulence is much less problem-dependent than the large 
scale turbulence so that the sub-grid-scale turbulence can 
be represented by relatively simple models. This approach 
appears very promising for solving three-dimensional time- 
dependent problems. Kim (1985) used 128x129x128 
computational grid points for direct numerical simulation of 
turbulent flow in a channel. He used the temperature as a 
passive scalar to investigate the interaction of the wall- 
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layer structure with the outer layer. Herring (1979) gives a 
short introduction and overview to subgrid modelling and 
Love and Leslie (1979) give more details. The greatest 
drawback with this approach is the huge amount of computing 
time involved. For this reason, they are being looked upon 
as methods of improving the modelling approximations of 
simpler closures. 

In view of the similarity of the density-weighted 
averaged equations to those for uniform-density flow, the 
idea has emerged that existing uniform density models may be 
adapted to non-uniform density turbulent flows simply by 
substituting the density-weighted variables for the standard 
average variables in a particular model. Bilger (1976, 1977, 
1979) advocates this hypothesis of turbulent model 
similarity. However, not much work has been done in this 
area. Therefore, our present objective is to try to 
establish the validity of constant-density models for 
variable density flows. 

1.2 OBJECTIVES 

The main objectives of this research are: 

A. To evaluate and identify the most general and efficient 
model for turbulent momentum exchange in swirling and non- 
swirling combustor flow calculations. The evaluations will 
include k-e model, algebraic stress models and full Reynolds 
stress models. 


13 


B. To provide a review of the existing turbulent scalar flux 
models and to evaluate and identify the most general and 
efficient scalar flux model for swirling and non-swirling 
combustor flows. 

C. The validity of these models for variable-density flows 
is examined and their applicability is demonstrated by 
comparison with measurements. 

D. The turbulence models identified above will be applied to 
calculate premixed and non-premixed reacting flows using 
both fast and finite rate chemistry models. Their validity 
will be examined in detail. 

E. Finally, the effect of heat release on the turbulent flow 
field is examined, and the validity and extent of constant- 
density turbulence models for reacting flow calculations is 
assessed. 

1.3 OUTLINE OF THE REPORT 

The remainder of the report consists of five chapters 
and the accompanying appendices. 

In chapter 2, the problem of calculating turbulent 
flows is posed more precisely by introducing and discussing 
the density-weighted averaged equations governing the mean- 
flow quantities. The appearance of turbulent transport terms 
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in these equations makes apparent the necessity of 
introducing turbulence models. The heart of this chapter is 
the actual review of models in the sections 2.3 to 2.4; the 
models are discussed in order of increasing complexity. 
Section 2.5 considers the extension of the turbulence models 
to include the low Reynolds number region that is always 
found in the immediate vicinity of a smooth wall. Although 
this region usually occupies less than 1% of the flow 
domain, it is of significant importance because 50% of the 
change in mean velocity occurs across it. Section 2.6 turns 
the attention to the modelling of combustion processes. It 
will discuss the cases of non-premixed and premixed 
reactants, respectively. 

Chapter 3 presents the details of the solution 
procedure adopted for the highly coupled and non-linear 
governing equations. This chapter briefly discusses 
numerical (false) diffusion which may or may not seriously 
affect the accuracy of the solution. A scheme is introduced 
to reduce this source of error. 

In chapter 4, the effects of four different pressure- 
strain models on three different algebraic stress closures 
(ASM) for swirling and non-swirling turbulent flows are 
investigated. The results are compared with the standard k-e 
model. Having demonstrated the effect of the pressure-strain 
correlation and determined its suitable model, attention is 
turned to the performance of the Reynolds stress closure 


15 


(RSM) . A comparison is made of the performance of a high and 
a low Reynolds number Reynolds stress closure. Effects of 
two different turbulent diffusion models on RSM are also 
analysed and the results are compared with the standard k-£ 
model and some algebraic stress closures. For swirling 
flows, comparison of standard k-e model, ASM and RSM are 
made using the diffusion model and pressure-strain model 
which are found to perform the best for non-swirling flow 
calculations. As for the scalar field calculations, effects 
of two pressure-scalar-gradient models on two different 
algebraic scalar flux models (AFM) are analysed and the 
results are compared with the full Reynolds stress/flux 
closure. Finally, comparison of three different ASM and the 
k-€ model are made for variable-density swirling flows and 
the results are discussed. 

Chapter 5 discusses the results of non-premixed and 
premixed combustion models. Both fast and finite rate 
chemistry models are applied to non-premixed combustion. 
However, in contrast to non-premixed flames, premixed flames 
require the consideration of finite rate reaction only. In 
this case, a two-step reaction process is used and the mean 
formation rates are calculated from the Arrhenius reaction 
rates and also the eddy break-up model. 

Finally, chapter 6 summarises the main conclusions 
emerged from this study and put forward some recommendations 
for future work. 
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CHAPTER 2 


GOVERNING EQUATIONS FOR VARIABLE-DENSITY FLOWS 
2.1 MEAN EQUATIONS IN FAVRE -AVERAGED FORM 

In this section, the equations which govern the 
distribution of the mean flow quantities are presented. 
These equations are derived from the conservation laws of 
mass, momentum and scalar which can be expressed in 
Cartesian tensor notation as: 

mass conservation : 


d£ 

dt 


+ 



1 


( 2 . 1 ) 


momentum conservation : 


du. , ~ <9u. 

■at' + - 


-SB. + 2t, j 
3x. dx! 


scalar conservation : 


dd ~ do 

P TT + pu. 0 — 
dt y J dx. 


= - ^ 


dx. 


( 2 . 2 ) 


(2.3) 
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where u. is the component of instantaneous velocity in the 
x. direction, p is the instantaneous static pressure, 6 is 
an instantaneous scalar quantity, is the diffusion flux 

of scalar 6 in the x. direction, and r fJ is the 
instantaneous viscous stress tensor. 

Turbulent flows contain motions which are much smaller 
than the extent of the flow domain. In order to resolve 
these motions in a numerical procedure, the mesh size of the 
numerical grid would have to be even smaller. Storing the 
flow variables at so many grid points is beyond the capacity 
of present computers. For this reason, turbulent flow is 
normally analysed using statistical methods. Following 
Reynolds, the instantaneous quantities are separated into 
mean values and their fluctuations. For variable density 
flows two types of decomposition can be used; either the 
unweighted form conventionally used for constant density 
flows or the density-weighted decomposition suggested by 
Favre (1969) . The unweighted decomposition and averaging of 
flow variables are represnted by 

ui = Ui + 

where 


U 


i 


= Lim | 

T — ► OO 


T+t 



t 


dt 


and u| = 0 
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are 


which the density-weighted decomposition and averaging 
defined as: 


Uj 


= U i + u A 


where 


T+t 


u i = Lim £ 




fiUf dt 
p 


oo 


Since by definition pu.=0, it can be shown that u.*0 and 
pu.=-p'u\. The averaging time T is long compared with the 
largest turbulence time scales, but shorter than the period 
over which the average flow quantities may vary. In general, 
the averaging process should involve ensemble averaging, but 
for stationary flows, time averaging and ensemble averaging 
are the same (Lumley, 1970) . For high Reynolds number flows, 
after noting that density weighting is not to be applied to 
either the pressure or density , the Favre-averaged forms of 
the equations of continuity and conservation of momentum and 
scalar may be written as 
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(2.4) 
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( 2 . 6 ) 


The unweighted averaging technique results in similiar 
equations but with two differences. Firstly, the density 
weighted quantities are replaced by unweighted quantities. 
Secondly, correlations such as p u and p d appear in the 
equations . 

There are some advantages for selecting the Favre 
averaged forms of the equations. It provides equations 
describing the mean values of those quantities which are 
conserved. For example, the component of mean momentum per 
unit volume in the x^direction is pU. and not pU,.. In 
addition, the equations which arise from the application of 
the density-weighted averaging technique are of a much 
simpler form and are more easily interpreted than those 
obtained from unweighted averaging. For the reasons stated 
above, we will concern ourselves only with the density- 
weighted averaged forms of the governing conservation 
equations from now on. 
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As a consequence of the non-linearity of the equations 
(2.2) and (2.3), the averaging process used introduces 
unknown correlations. As a result, the equations are no 
longer closed and closure assumptions are necessary before 
solution is possible. The problem of calculating these 
correlations can be approached at different levels of 
complexity. It is the intention of this chapter to discuss 
in detail the appropriate turbulence closures for variable 
density flows. These include both turbulent viscosity and 
full Reynolds stress models. 

2.2 REYNOLDS EQUATIONS IN FAVRE-AVERAGED FORM 

For high Reynolds number turbulent flows, the equations 
that govern the transport of the Reynolds-stress tensor may 
be concisely expressed in Cartesian tensor form as 


p |t< u i u j> + exl u i u i ) = 
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-ax k u i u j u k) 
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(2.7) 
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Equation (2.7) is obtained by subtracting equation (2.5) 
from (2.2), multiplying the result by u^ and adding to it 
the same equation with subscripts i and j interchanged, then 
averaging the resultant equation. It expresses the rate of 
change of the Reynolds stresses as the sum of four different 
processes: (a) diffusion by the action of velocity 
fluctuations, (b) transport by fluctuating velocity-pressure 
gradient correlations, (c) generation of Reynolds stresses 
through the interaction of stresses with the mean strain 
rate, and (d) dissipation of Reynolds stresses through 
viscous action. The contraction of this equation, that is 
when the equations for the three normal stresses (i=j=l,2,3) 
are summed up, yields the turbulent kinetic energy equation. 

The density-weighted averaged form of the equation for 
the scalar flux, may be obtained by multiplying the 
fluctuating velocity equation by the fluctuating scalar, 0; 
then adding to it the equation for 6 multiplied by u. and 
averaging the result. Neglecting the influence of buoyant 
force and viscous diffusion, the result is expressed as: 
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P § t (U,.*) + P U m ^(V)- -g^Cp UiV + p'tf 5 im ) + p'f x . 
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(B) 


-/ — — 9U. ^ ae x 

-p ( u^O — 1 + u- u w ' — ) 
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ij ax j 


- j - ^ ) 
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( 2 . 8 ) 


(C) 


(D) 


The physical interpretion of the terms on the right hand 
side of (2.8), from left to right, is: (A) the diffusive 

transport due to velocity and pressure fluctuations, (B) 
transport due to the fluctuating pressure-scalar-gradient 
correlations, (C) the generation of scalar flux due to the 
mean gradients of the velocity and scalar fields, and (D) 
viscous dissipation. 


2.3 MODELLING OF THE REYNOLDS EQUATIONS 

2.3.1 MODELLING OF THE u. u -EQUATION 

Of the terms on the right-hand side of equation (2.7), 
only the production term (c) can be calculated directly in a 
flow field where the velocities, Reynolds stresses and 
dissipation rate are known. Others, however, need to be 
modelled because they either include higher order 
correlations (term a) or correlations between turbulence 
quatities that are not known. The objective of this section 
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is to discuss each of these terms separately and to suggest 
approximate models for closure. In order to accomplish this 
task we shall make extensive use of models developed for 
uniform density flows and assume that the use of density- 
weighted averaging will partially account for the influence 
of density fluctuations (Jones, 1980) . 

The viscous dissipation correlation (d) is the least 
difficult term to model. This term represents the 
destruction of Reynolds stresses through viscous action. At 
high Reynolds number, the small scale structure of 
turbulence tends to be independent of any orientation 
effects introduced by the mean shear. Under these 
conditions, the small scale strucure of turbulence is in a 
state of local isotropy and the dissipation correlation as 
suggested by Kolmogorov (1941) can be expressed as 



' du. ' Su. 
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ki dx k kj Sx k 
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is the dissipation rate of turbulent kinetic energy per unit 
volume. This model introduces another unknown into the 
problem; namely e, the turbulent kinetic energy dissipation 
rate. Therefore, an equation for the dissipation rate must 
be solved in conjunction with the set of Reynolds-stress 
equations . 

Term (b) is conventionally partition into two parts 


u j 


dxj 


+ Uj 


Si ax k (5 jk p ' u i + p ' u j ) 


- P'( 


au. 


. . ^j 
axj ax^ 


i + 


( 2 . 10 ) 


where the first term on the right hand side of (2.10) is 
conventionaly identified as diffusive transport of uTUj due 
to p' and the second term is usually referred to as the 
pressure-strain term. This term is redistributive if the 
flow is incompressible, because the trace of this term is 
identically zero. 

For variable density flows, the conventional 
partitioning of (2.10) does not result in a redistributive 
term since aUj/ax^O. An alternate partition which gives a 
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redistributive term has been suggested by Lumley (1975 b).It 
can be written as 
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( 2 . 11 ) 
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It can be seen that the first term on the right-hand-side of 
(2.11) gives a zero trace whether the flow is incompressible 
or not. The primary function of this term is to change the 
relative levels of the normal stresses, and to redistribute 
the turbulence energy in the stress equations. 

From a practical point of view, all these arguments are 
quite immaterial, since none of the terms can be measured 
(Jones 1980) . However, from the view point of turbulence 
modelling, the difference may be important. Thus if we want 
to apply the constant density models (which are 
redistributive) for variable density flows the latter 
decomposition (2.11) is recommended. 

In order to gain further insight into the modelling of 
the redistributive term, a Poisson equation for the 
fluctuating pressure is derived. This is done by taking the 
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divergence of the equation for the fluctuating velocity u < . 
The result is expressed as: 


2 2 

d p • _ r 
dxt ” at 2 
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(pUj^U j —PUj^U j +pu i Uj+pUjU i ) 
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ax^axj 


( 2 . 12 ) 


Equation (2.12) is much more complicated than its constant 
density version because aUj/ax^O. According to Jones (1980) 
some of the influence of density fluctuations could be 
partially accounted for by considering the Favre-averaged 
equations rather than the Reynolds-averaged equations. 
Besides, since the complete modelling of variable-density 
flows is, at present, hindered by a lack of credible data, 
therefore, it is more practical to apply constant density 
models. 

Following Chou (1945) the redistributive part of the 
velocity-pressure gradient correlation, n^j , is divided 
into two types of terms: one containing fluctuating 
quantities alone (IIjj n ) and the second due to mean strain 
effects, (n^j ->) . Each contribution is ususally modelled 
separately. The term IT^j has long been recognized as the 
only mechanism in the stress transport equation (2.7) that 

I 

could promote a return to isotropy. Rotta s (1951) proposal 
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for the modeling of has generally been adopted and is 

used in most Reynolds-stress equation models: 
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ij / I 


— c, p 
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1 P k (u i u j " 25 ij k/3) 


(2.13) 


where e/k defines a time scale and is a constant to be 
determined from experimental considerations. 

A model to represent the interaction of mean strain and 
fluctuating velocities was first devised by Rotta 
(1972). In the absence of external forces this can be 
written as: 


n i j 1 2 " C 2 P k ( s ij " | 5 ij f^") 
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Later, Naot et al. (1973) proposed that the effect of mean 
strain rate is to promote a similar return to isotropy in 
the Reynolds stress production tensor, i.e. 
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where 


p ij s - * < u i u * gj, + u j“k ^ ) - 


p k . 2ii. 


Launder, Reece and Rodi (1975) suggested a model which 
into account effects due to both the symmetric 
antisymmetric parts of the rate of strain tensor, 
model takes the form 
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The remaining part of the fluctuating velocity-pressure 
gradient correlations is not redistributive and can be 
written as 
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m ax 
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(2 . 17) 


Bilger (1976) has argued that the first term is associated 
with noise generation and is small. In constant density 
flows and in situations where the Boussinesq approximation 
is appropriate this term vanishes. The second term can be 
lumped into the diffusive transport terms. Hence, the 
diffusive transport of turbulent stresses includes three 
mechanism: (1) transport through fluctuating velocities, (2) 

transport through pressure fluctuations, and (3) molecular 

transport, which is negligible at high Reynolds number. 

. ' ... 
Irwin s (1973) measurements of self preserving wall jets in 

an adverse pressure gradient show that the transport due to 

pressure fluctuations is negligible compared with that due 

to velocity fluctuations. Within these assumptions, the 

diffusive transport term reduces to the triple-velocity 

correlation only. Daly and Harlow (1970) modelled this term 

as 
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(2.18) 


' u i u j u k 


= C - Ui,u_ — (uu- ) 

S € k u m dx m ' ' J J 


The left-hand-side of equation (2.18) is tensor-invariant 
but the right-hand-side alters under permutation of the 
indicies i,j and k. In order to have a model that is 
consistent tensorially, Hanjalic and Launder (1972) propose 
the following form for this correlation 


-UiUjUk 


C s"J u k u m + u i u m ax m (u jV + u j u m dx m < u i u k> ] 

(2.19) 


m 


m 


This formulation contains the Daly and Harlow modelling as 
its first term and possesses the correct tensor properties. 
Both (2.18) and (2.19) imply gradient transport of Reynolds 
stresses but whereas in (2.18) the dif fusivities are the 
same for each stress component, the dif fusivities evaluated 
from (2.19) are different for each stress component. 


2.3.2 MODELLING OF THE u 9 -EQUATION 

The mean field production term (C) in equation (2.8) is 
due to the combined action of mean velocity and mean scalar 
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gradients, the former tending to increase the scalar 
fluctuations and the latter the velocity fluctuations. This 
term is exact and does not need to be modeled. At high 
Reynolds number, the fine scale motions are assumed to be 
isotropic. Since there is no first-order isotropic tensor, 
the dissipative term (D) is zero. Consequently, model 
approximations have to be proposed for the diffusion term 
(A) and for the fluctuating pressure-scalar-gradient 
correlation term (B) . 

By neglecting pressure transport, Launder (1976) 
proposed the following gradient approximation for term (A) : 


- u iV = C S0 7 [u m u k ax + u i u k a'x k (u ^ )] 


( 2 . 20 ) 


where the coefficient C s q has been taken as 0.11 in 
conformity with the value of C s in equation (2.19). 

The most important process requiring approximation is 
the pressure-scalar-gradient correlation n^. This term is 
the counterpart of the pressure-strain term in the Reynolds 
stress transport equations. This is the main term to 
counteract the production of the scalar flux and thus to 
limit its growth. With the aid of the Poisson equation for 
the fluctuation pressure, two different mechanisms which 
contribute to can be identified. One is due to pure 
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turbulence interaction, and the other arising from mean 
strain effect. Therefore, can again be partitioned into 


n i0 n i*,l + n i0,2 


( 2 . 21 ) 


The most widely used model for the first part is due to 
Monin (1965) 
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i0, 1 



( 2 . 22 ) 


Equation (2.22) is a direct counterpart of Rotta ' s return to 
isotropy approximation (2.13) for the pressure-strain term 
^ij,l* Analogous to the pressure-strain model (2.15), 
Launder (1975) proposed that the second part should be 
proportional to the scalar flux production due to the mean 
strain, or 
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(2.23) 


2.3.3 THE DISSIPATION RATE TRANSPORT EQUATION 

In anticipation of the fact that the quantity e is 
required in turbulent flow calculations, we need to solve a 
transport equation for e. For high-Reynolds-number variable 
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density flows, the transport equation for turbulent kinetic 
energy dissipation rate is given by Jones (1980) as: 
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Additional terms could be added to equation (2.24) for 
variable density flows. These terms are: 


e £lu ; ap 1 
k p ax^ 


- 9U- 
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However, this introduces more empiricism into the equation 
and there are no reliable means to determine the additional 
empirical constants. As a result, the additional terms for 
(2.24) are commonly neglected. 

Equation (2.24) shows that the rate of change of e is 
balanced by the diffusive transport, the generation of 
vorticity due to vortex stretching connected with the energy 
cascade, and viscous destruction of vorticity. 
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2.3.4 THE SCALAR FLUCTUATION TRANSPORT EQUATION 

In connection with variable-density flow calculations, 
it is necessary to determine the scalar fluctuation 
quantities (Spalding, 1971) . For high-Reynolds-numbers , the 
transport equation for the scalar fluctuation is given by 
Jones (1980) as: 
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The equation expresses the fact that the level of ~¥ 
following a mean stream-line will change through an 
imbalance of the generation rate of the scalar fluctuation 
by gradient in © (AA) , through diffusive transport produced 
by turbulent velocity fluctuations (BB) and the dissipation 
of fluctuations due to molecular diffusion in the fine scale 
motions (CC) . Equation (2.25) is the most straightforward of 
all turbulence transport equations; it resembles the 
turbulence kinetic energy equation except that pressure 
transport is absent from (2.25). 
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2.4 DIFFERENT LEVELS OF CLOSURE MODELS 

The task of turbulence closure is to determine the 
turbulent flux terms in the mean-flow equations so that 
these can be substituted into the equations and thus allow 
the velocity, temperature and concentration fields to be 
solved. Turbulence models may be classified according to 
whether the turbulent momentum and heat/mass fluxes are 
assumed to be locally determined or whether they are 
obtained from their transport equations. The former approach 
is known as mean flow closure while the latter is called 
Reynolds stress closure. 

In mean flow closures, turbulent viscosity relates the 
Reynolds stresses to the local strain field by an expression 
similar to that for laminar flow. However, the turbulent 
viscosity is a function of the properties of the turbulence 
and not of the fluid. A wide variety of models are available 
in the literature. However, these have only been applied to 
particular classes of thin shear flows. Since these models 
are not discussed in detail in later sections, it may be 
helpful to briefly mention some of the more pertinent work 
here. The simplest model that provides information about the 
distribution of the eddy viscosity, is the mixing-length 
hypothesis, generally attributed to Prandtl (1925). Prandtl, 
stimulated by kinetic gas theory, assumed that the eddy 
viscosity is proportional to a mean fluctuating velocity 
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t? and a mixing length £ m . Considering shear layers with only 
one significant turbulent stress (uv) and velocity gradient 
(dU/dy) , he then postulated that t? is equal to the mean 
velocity gradient times the mixing length. The closure 
problem is then reduced to a determination of the mixing 
length l m . Prandtl found that simple algebraic prescription 
l m leads to fairly good results for simple shear flows. 

Recognizing that the level of turbulence at a point in 
the flow is determined not just by local events but by 
convective and diffusive behavior too, Prandtl (1945) 
proposed that the square root of the turbulence kinetic 
energy, k should be adopted as the characteistic fluctuating 
velocity and that its value should be determined from a 
transport equation. The mixing-length £ m is typically 
described with the aid of empirical information. The general 
form of the equation he proposed is still in use today in 
models of the "one-equation" type. Two further one-equation 
models also deserve mention. Nee and Kovasznay (1969) 
devised a transport equation for the turbulent viscosity 
itself while Bradshaw et al. (1967) developed a simplified 
transport equation for the shear stress. The latter model 
was suitable for calculating thin shear flows where the 
shear stress did not change sign; it has been extended to 
allow reversal in the shear stress across the layer 
(Bradshaw et al. , 1974). Unlike boundary-layer flows on a 
smooth wall, there are no simple rules that can serve to 
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give the length scale in a separated flow. The only way to 
obtain a reasonably accurate estimate of the length scale 
distribution is to calculate it, like the turbulent kinetic 
energy, from a simulated length scale equation. 

Two-equation models are the simplest models that are 
suitable for complex flows in which the length-scale cannot 
be prescribed empirically in an easy way. The various 
length-scale equations proposed in the literature have very 
similar performance? the e-equation become most popular 
because of its relative simplicity. The standard k-e model 
is based on the assumption that the eddy viscosity is the 
same for all Reynolds stresses ujuj . However, in certain 
flow situations the assumption of an isotropic eddy 
viscosity is too crude; for example it does not produce the 
turbulence-driven secondary motions in square ducts which 
have been observed experimentally (Meroney, 1976) . To allow 
for the non-isotropic behavior of the eddy viscosity and to 
account for the effect of body forces (buoyancy , rotation) , 
the k-e model is refined by introducing algebraic stress 
model. This model is based on a simplification of the 
Reynolds stress equations that allows the equations to be 
reduced to algebraic equations for u^j . Since the 
quantities k and e are present in these equations, their 
transport equations also have to be solved. 

In the Reynolds stress closure, the stress components 
are no longer coupled to the local mean strain rates but are 
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determined from their transport equations. Models employing 
transport equations for the individual turbulent stress and 
flux components simulate the turbulent processes more 
realistically and are therefore potentially more general 
than the simple models. 

Combustor flows are very complicated. They contain 
different flow regions such as: (1) inviscid core region, 
(2) recirculation zone, (3) shear region between regions (1) 
and (2) , (4) reattachment region, (5) recovery region and 
ordinary boundary layer along the wall. Each region has its 
particular turbulence characteristics. Therefore, any 
proposed model has to be general enough to predict all these 
different regions if it is to be successful in combustor 
flow calculations. It is not possible to specify empirically 
the length-scale distribution in such complex flows. As a 
result, attention is given to turbulent transport models in 
which at least one turbulent velocity and length-scale is 
found from the solution of approximated transport equations. 
Only turbulence models of at least this degree of complexity 
are suitable for complex combustor flow calculations. 
Therefore, the standard k-e model, algebraic stress model 
and Reynolds stress model are the only qualified models for 
combustor flow calculations. 
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2.4.1 k-e MODEL 


The k-e model developed by Jones and Launder (1972) is 
the simplest model which is suitable for recirculating flow 
calculations. This model achieves closure by using a 
gradient transport model for uTCt , or 
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The effective (turbulent or eddy) viscosity n t 
appearing in (2.26), is defined in terms of a characteristic 


length and 

velocity. If 
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approximated by k 1/2 , fi t can be expressed as 



where C^is a constant of proportionality. 

The k-e model is the simplest model that allows the 
characteristic length scale of a wide range of complex flow 


fields to 

be 

determined. It has 

been 

applied 

to 

the 

calculation 

of 

various free shear 

flows 

(Launder 

et 

al„ 


1973) , recirculating flows (Sindir, 1983a, 1983b) , and 
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confined reacting flows with and without swirl (Smith and 
Smoot, 1981) . In direct comparison of the k-e model 
predictions with experimental data, the model is successful 
in predicting the basic features of turbulent flows. When 
significant streamline curvatures are introduced into the 
flow field, the k-e model does not adequately account for 
the enhanced turbulence diffusion caused by the extra strain 
rates associated with streamline curvature. The standard k-e 
model is based on the isotropic eddy viscosity assumption. 
However, in the calculation of simple thin shear layers 
where only the shear stress, uv, appears in the equations, 
the closure is not influenced by the isotropic assumption. 
In other flow situations where such simplifications are not 
realized, models based on isotropic turbulent viscosity 
often lead to unsatisfactory flow predictions. To allow for 
the non-isotropic behavior of the eddy viscosity in such 
cases, the k-e model is often modified by introducing the 
algebraic stress model to replace the isotropic eddy 
viscosity assumption. The axisymmetric form of the turbulent 
flow equations is given in appendix A for the k-e model. 

2.4.2 THE ALGEBRAIC STRESS/FLUX MODELS 

The algebraic stress model (ASM) is a special case of 
the Reynolds stress model. It relates the individual 
stresses to mean velocity gradient, turbulent kinetic energy 
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and its dissipation rate by way of algebraic expressions. 
Algebraic stress models, hereafter denoted by ASM for short, 
can be classified into two categories. The first is based on 
a local equilibrium assumption for the turbulence field, 
whereby the turbulence transport terms are neglected 
compared to the local production and dissipation of 
turbulence. A second class of ASM is based on the local non- 
equilibrium assumption. Approaches of this kind, where the 
convection and diffusion transport of turbulent stresses are 
approximated, have been developed by Mellor and Yamada 
(1974,1982) and Rodi (1976). 

The equilibrium ASM, hereafter denoted by ASM/E, is 
obtained when the convective and diffusive transport terms 
in the Reynolds stress transport equation (2.7) are assumed 
to be zero. With this simplification, (2.7) reduces to 



| S . . e +n . . 
3 ij lj 
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With suitable modelling for n^j , (2.28) becomes algebraic in 

u,-u. and can be solved in terms of k, e and S fj . . Many flows, 
in particular thin shear layers, evolve so slowly that 
neglecting the convection and diffusion of turbulent 
stresses and fluxes can be justified. It should be mentioned 
however, that these simplified relations lead to an 
inconsistancy in the normal stresses when dissipation is not 
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balanced by production, because the resulting normal 
stresses do not sum up to 2k as they should. 

Mellor and Yamada propose to simplify the Reynolds 
transport equations through an order-of-magnitude argument 
based on a 2 = Ofa^) 2 , where a^. is a non-dimensional measure 
of anisotropy and is given by 
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The order-of-magnitude argument is performed on an equation 
for (2a { .k) which is obtained by subtracting the product of 
Sjj/3 and the transport equation for turbulent kinetic 
energy from the transport equation for uTCTj. . The resultant 
equation becomes 
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where the differential operator S£ is used to denote the 
combined convective and diffusive transport operators. Terms 
in (2.30) are evaluated in powers of a and terms of order a 2 
and higher are neglected. The result is 
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Rodi's (1976) approach is to eliminate the transport 
terms in the governing equations by assuming 


# (UjUj) 


(HiHj. ) #(k) 
k 


(2 . 32) 


This assumption is correct when uTuj/k is constant over the 
whole flow field. However, its validity for flows with 
varying uTUj/k ratios has not been demonstrated. Another 
point to note is that u]Uj/k is constant only for flows that 
are in local or near equilibrium (Mellor and Yamada, 1982). 
Nevertheless, this assumption leads to 

^j( p k - PO = KPij - §S i;j ie + n i:J ) (2.33) 


which gives a set of algebraic equations for u.Uj after n^j 
has been appropriately modeled. 

The approaches of Mellor and Yamada, and Rodi are 
essentially similiar to the same order a 2 . This can be seen 
by considering assumption (2.32), which can be written with 
the help of (2.29) as 

*<*4^ " | 5 ij k > = #(*) k) (2.34) 

k 
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Since ^(k) = P k - pe , and by virtue of (2.29) 


& ( u i u j " | 6 ij k ) = 2 ( p k " P Oa i j (2.35) 

The term 2(P k -pe)a jj . is of order a 2 and (2.35) would again 
lead to the results of (2.31). Rodi's assumption also leads 
to the neglect of terms of order a 2 in the Reynolds 
transport equations, and the two different approaches are 
essentially similiar to the same order, a 2 . 

If the two non-equilibrium ASM are denoted by ASM/MY 
and ASM/R, then ASM/MY is given by (2.31) and ASM/R is given 
by (2.33). The three ASM yield a set of algebraic equations 
for u i Uj once a suitable model is proposed for the pressure- 
strain terms, IT^j . Therefore, the performance of the ASM 
will be influenced greatly by the model adopted for n^j . 

A similar approach can be applied to model the scalar- 
flux equations. As a result, no explicit statement about the 
form of diffusive transport in the scalar flux equation 
(2.8) is needed. The resultant model is the scalar-flux 
counterpart of the algebraic-stress-modelling (ASM) method. 
According to the order-of-magnitude analysis of Mellor and 
Yamada (1974) the turbulent scalar equation (2.8) reduces to 
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where local isotropy of fine scale turbulence have been 
assumed. Note that (2.36) does not imply local balance of 
production and dissipation of 6 2 . 

Gibson and Launder (1976) proposed to approximate the 
turbulent scalar flux equation in a parallel manner to 
(2.33). Their result gives 

<p « -*» > 

(2.37) 

where R and e. are respectively the production and 

U U 

dissipation of the scalar fluctuation 6 2 . Gibson and Launder 
(1976) then argue that since there is only a weak coupling 
between the d 2 and u70 equations the assumption of local 
equilibrium for 6 2 is suffeciently accurate. Then (2.37) may 
be further simplified and written as: 


"tt 2Ul* 


+ 3 X (P u iV + s im p # > - 


m 


ui . u.0 

— 1 - (Pk-PO - — ^7 
2k /)2 


- . 9© 

(u i u j 


— 7- 9U. 

Uj.6 — 1 

3 9X4 


) + P 


86 

dXi 


+ 


~ U; 
P — L 


, 9U. , , 

2k lUj * X j ° 


= 0 


(2 . 38) 


46 


Relations (2.36,2.38) and (2.31,2.33) together with the k- 
and e — equations form a two— equation model which can account 
for the different production or destruction processes acting 
on the Reynolds stress and flux components, and to a certain 
approximation also for the transport of these components. 
This model is therefore cosiderably more general than the 
standard k — e model which employs an isotropic eddy 
viscosity. These governing equations in the axisymmetric 
coordinate system are given in appendices B through D. 


2.4.3 REYNOLDS STRESS MODELS 

The models reviewed so far assume that the local state 
of turbulence can be characterized by one velocity scale and 
that the individual Reynolds stresses can be related to this 
scale by the eddy-viscosity expression. This relation often 
implies that the transport of the individual stresses is not 
adequately accounted for, even if the transport of the 
characterizing velocity scale is. In order to allow for the 
different development of the various Reynolds stresses 
representing various velocity scales in complex flows and to 
account properly for their transport, models which employ 
transport equations for the individual stresses uTu. must be 
applied. Analogous transport equations are required for the 
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turbulent heat/mass fluxes uTfl, and the models based on 
these equations are often referred to as Reynolds stress 
models. The derivation of u^Uj" and equations (2. 7, 2. 8) 

and their modellinq to obtain a closed system were explained 
in previous sections. In general, there are six components 
of the Reynolds stresses and three components of the 
turbulent scalar fluxes. The modelled equations for these 
components, the turbulence dissipation rate e, the scalar 
fluctuations 6 2 and the mean flow equations (2. 4-2. 6) allow 
the calculation of non-reacting variable density flows. 
However, for turbulent reacting flows, a combustion model 
must also be added. A complete set of the Reynolds 
equations in cylindrical coordinate are presented in 
appendices E and F. 

2.5 NEAR-WALL FLOW MODELLING 
2.5.1 WALL FUNCTION 


At a solid boundary the no-slip condition applies. This 
means that both mean and fluctuationg velocities are zero. 
In contrast, the dissipation rate e is finite and requires 
special attention. When the boundary conditions are 
specified right at the wall, the equations must be 
integrated through the viscous sublayer. This is undesirable 
because the high-Reynolds-number turbulence models 
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introduced above are not applicable in this region. However, 
integration through the sublayer can be avoided by using the 
law-of-the-wall function to connect the wall conditions to 
the dependent variables just outside the viscous sublayer. 
The value of k at the first grid point which is selected to 
be outside the viscous sublayer is determined from the k- 
balance for a control volume adjacent to the wall and with 
diffusive and convective transport through the wall set 
equal to zero. On the other hand, e at the first grid point 
is calculated from the equilibrium condition. When turbulent 
transport equations are solved, the individual stresses 
(relative to friction velocity) are specified at the first 
grid point (Launder et al. ,1975). More discussion regarding 
the wall functions are presented in appendix G. 


2.5.2 DIRECT CALCULATION 

All the models described so far imply negligible 
viscosity effect on the energy containing motions and 
negligible effect of the mean strain field on the 
dissipative phenomena. These assumptions are generally valid 
for high Reynolds-number turbulent flows except very close 
to a wall. Here, the no slip condition at the wall ensures 
that viscous effect will always be important in the 
immediate vicinity of the wall. Although the viscosity 
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affected zone has a thickness of two or more order of 
magnitude less than the characteristic dimension of the 
flow, almost 50% of the velocity change from the wall to the 
free stream occurs in this region (Launder , 1984) . 

Fortunately, it is practical to ignore many of the 
complexities inside the viscosity dominated layer by noting 
the fact that the important mean and turbulent flow 
quantities highly depend on the normal distance from the 
wall. In high-Reynolds-number modelling, the calculations 
are carried out to the vicinity of the wall and all the 
dependent variables are then matched to their corresponding 
values determined from approximate wall functions. 

Transpiration through the wall, swirl, steep temperature or 
streamwise pressure gradients are just a few examples that 
may cause the near-wall region to differ from its so called 
universal behavior. To account for these various influences 
it is necessary to extend the calculation right up to the 
wall itself. 

In order to provide predictions of the flow within the 
viscous sublayer, the form of the models given for high 
Reynolds number flows must be modified in three ways. These 
are: (1) retain the viscous diffusion of k, e and uTuj, (2) 

constants in e and eddy viscosity equations become dependent 
on the turbulence Reynolds number, (R T =k 2 /i/€) , (3) further 

terms are required due to the fact that the dissipation 
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processes are not isotropic and this effect must be 
accounted for. 

Several two-equation, low-Reynolds-number turbulence 
models have been proposed by different researches (Launder 
and Sharma, 1974; Hassid and Poreh, 1978; Chien, 1980; 
Reynolds, 1976; Lam and Bremhorst, 1981). A systematic 
evaluation of these models has been performed by Patel et 
al. (1985). From an overall examination of the results for 
all the test cases, the model of Chien (1980) performs 
relatively better than the others. Therefore, this model has 
been selected here to calculate the turbulence energy and 
dissipation rate. The final modified form of the k and e 
equations proposed by Chein (1980) are 
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where 


f 2 = 1 - | e" R T/3 6 


and y is the normal distance from the wall. 
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Recently, So and Yoo (1986) applied asymptotic analysis 
technique in a way parallel to Chien (1980) to modify the 
viscous dissipation term in the Reynolds stress transport 
equations to account for viscous behavior near a wall. 

2.6 TURBULENT COMBUSTION MODELS 

Turbulent reactive flows are difficult to predict for 
the simple reason that neither turbulent transport nor 
chemical kinetics are adequately understood. Specific 
problems of particular interest are the evaluation of the 
mean formation rates and the effects of heat release on the 
turbulence structure and, in turn, on the turbulence flux 
closure models. In order to illustrate the effect of 
turbulence on the mean formation rates, consider a simple 
irreversible reaction such as 

kl 

A + B — C (2.41) 

The instantaneous formation rate of C is represented by the 
Arrhenious equation, or 

R C = k Q p 2 m A m B exp(-El/RT) (2.42) 
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where m^ is the instantaneous mass fraction of species i, El 
is the activation energy and k Q is the pre-exponential rate 
constant 

In order to evaluate the mean formation rate R c , it is 
necessary to decompose each of the fluctuating variables 
into their mean and fluctuating components. It is emphasized 
that the time-mean rate expression is not simply the 
reaction rate calculated from the time-mean variables. For 
example 

- _2 

R n = m m p2k exp(-El/RT) # m m p k exp(-El/RT) (2.43) 

Instead, after applying Reynolds decomposition and 

appropriate simplification we obtain: 


R, 


= p 2 m A m^ 0 exp(- ^ [i+ « + 2 S + 2^ + 


RT 
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m A m B 


P m A p m B 


El (aEJ + ( =1. + .. ] 


(2.44) 


RT m A T 


2RT 


T 2 


which is an infinite series in terms of m A m B , p'm A , m A T' , 
etc. For K c to be evaluated correctly, an infinite number of 
these moments containing the variables T, m A and m' B would 
have to be determined. This rate expression is for the 
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simple example given by equation (2.44). For any realistic 
combustion scheme, many of these equations have to be 
considered. It is ideal to be able to handle reacting flows 
in which combustion may occur via a large number of finite- 
rate reaction steps. However, this would necessitate the 
solution of conservation equations for the mean value of 
each of the independent species, which in turn requires the 
evaluation of the mean formation rate of each species. It is 
clear that closure of turbulent reacting flow is a 
formidable problem. 

At present, only PDF transport equation formulation 
offers the possibility of handling large numbers of reacting 
species (Pope, 1976), but in view of computer storage 
requirements, run times and multidimensionality of the 
approach, their suitability requires further investigations. 
Fortunately, most of the reactions associated with the high 
temperature oxidation of hydrocarbon and hydrogen fuels 
usually have time scales very short compared with those 
characteristics of the turbulence field and the assumption 
of fast chemistry thus provides a reasonable description of 
the equilibrium composition, temperature and mixture 
density. However, the estimation of unburnt fuel, formation 
of pollutants, and the study of such phenomena as ignition 
and blow out require consideration of finite rate reactions. 
The fast chemistry assumption can only be applied to non- 
premixed combustion systems where fuel and oxidant are 
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injected separately. However, in premixed combustion system, 
where fuel and oxidant have been mixed prior to reaction, 
only finite-rate chemistry assumption can be invoked because 
the species are allowed to come into contact without burning 
occurring at the same time. In view of these differences the 
following discussion of combustion models will be divided 
into non-premixed and premixed combustion. Under non- 
premixed combustion, both the fast and finite-rate chemistry 
models are discussed. 


2.6.1 FAST CHEMISTRY MODEL FOR NON-PREMIXED COMBUSTION 

A practical idealized approach in non-premixed 
combustion (two-feed) system is based on the assumption that 
the chemistry is fast. Its first principal feature is the 
neglect of all intermediate reactions so that pure fuel and 
pure oxidant will react to form the products the moment they 
are in contact. Its second main feature is the assumption 
that the effective diffusivity coefficient of all species is 
the same. The first is justifiable only by the need to 
simplify. The second cannot be very far from reality for a 
turbulent flow. 

When the chemistry can be simplified to a one-step 
reaction such as 
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1 Kg fuel + i Kg oxidant -► (1+i) Kg product 


we have 


R fu = Rox/i 

where R fu and R ox stand for the rates of creation by 
chemical reaction of fuel and oxidant Respectively. 

The consequence of these assumptions is that a suitably 
chosen linear combination of the equations describing the 
conservation of two-unpremixed reactants yields an equation 
whose form is identical to that describing the convection 
and diffusion of chemically inert species. The equation 
which results will have no source term. The dependent 
variable of the resultant equation is usually taken to be 
the mixture fraction and is defined as 

e = h ~ ^n) (2.45) 

(t'f 

where 7 is given by 

7 = ^fu — ^ox/^ - (2.46) 

or 

7 s m fu + m pr /i+l (2.47) 
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or 


7 = m ox + i m pr /l+i 


(2.48) 


Here itif u , m ox , it»p r , and i are the mass fraction of the fuel, 
oxidant, product and the stoichiometric oxidant requirement 
per unit mass of fuel. Subscript F and 0 stand for fuel and 
oxidant stream, respectively. 

Let us now presume that the reaction is every where 
complete. This gives the following conditions for 0: 


0 < © < © st 


m fu = 0 

m OX = ^OX,© ” © st ) 


(2.49) 


and 


© st < © < 1 


m 


fu 


= m 


m 


ox 


= 0 


(© -0,,) 

fu,F (1 -0^) 


(2.50) 


where © st is the stoichiometric value of the mixture 
fraction and is characterized by: 

@st = 1 + 1 ( m fu,F/ m ox,o) (2 ' 51) 

The mass fraction of the products of combustion can then be 
obtained from 
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fu 


- m 


ox 


(2.52) 


Equations (2.49,2.50) are the result of the fact that 
thermodynamic equilibrium is supposed to prevail throughout. 
This means that finite values of both fuel and oxidant 
concentration cannot prevail at the same point. 

For adiabatic operation of gaseous flame, the standard 
enthalpy is a conserved scalar and thus with the assumption 
of equal diffusivity of heat and mass, the local enthalpy, h, 
may be calculated directly from 


h = © h F + (1-©) h 0 


(2.53) 


where 
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T 



dT + H 
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m 


fu 


(2.54) 


H c and Cp are the heat of combustion and specific heat of 
the mixture at constant pressure, respectively. If the 
specific heat is assumed to be independent of temperature, 
then the temperature can be calculated as 


T 


h H c m fu 


(2.55) 
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As has already been discussed throughout this chapter, 
solution of the turbulent, variable-density, differential 
equations requires a properly averaged density. It can be 
seen from equations (2.4-2. 6) that the conventional mean 
density is required. However, according to these equations 
the Favred-averaging mixture fraction (scalar) is predicted 
along with the variance for the Favre-averaged probability 
density function (PDF) . By using the local instantaneous 
equilibrium approximation and convoluting over the PDF, the 
Favre-averaged density would result. On the other hand, if 
it is realized that the Favre-averaged of the inverse 
density is equal to the inverse of the Reynolds-averaged of 
the density, Then, 

1 = [ 1 £Xa & (2.56) 

p p 

0 

r-' 

where F(0) is the Favre probability density function of the 
mixture fraction and is unknown. A straightforward approach 
to determine this PDF is to assume a PDF with two 
parameters. The two parameters are determined from the first 
two moments of the PDF which are obtained by solving their 
respective transport equations. Different two-parameter 
PDF 1 s have been proposed and tested by a number of 
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researchers (Spalding, 1971 ; Lockwood and Naguib, 1975 ; 

Rhodes, 1975) . However, evidence in support of the p- 
distribution has been provided by Jones (1977) in his 
calculation of diffusion flames. We shall assume a ,0- 
distribution function for F(?) so that 
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where 0 < 6 < 1 and a > 0 , b > 0. 
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constants a, b can be found from (2.58) and (2.59) and are 
given by 


a = 8 [ 9 I 1 - SI 
0% 


1 ] 


( 2 . 60 ) 


= a — (_1 - ©) 
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(2.61) 
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Therefore, once 0 and 0 2 are known from the solution of 
their transport equations, a and b can be determined and 
consequently the mean density can be calculated at each node 
point. 

2.6.2 FINITE-RATE CHEMISTRY MODEL FOR NON-PREMIXED 
COMBUSTION 

An important problem in finite-rate chemistry is 
choosing an appropriate level of complexity, in view of the 
large number of species and chemical reactions taking place. 
One solution to this problem is the use of a global approach 
that reduces chemistry to the specification of an overal 
global oxidation scheme. Both m fu and m QX can have non-zero 
values at the same point at the same time. This model can 
predict quantities of interest: fuel consumption and heat 
release rates. However, it is the evaluation of the mean 
formation rate which presents problems. If fluctuation terms 
are neglected and the rate evaluated in terms of the mean 
values of the quantities involve in the reaction rate, the 
result may be in serious error. To reduce this error, the 
eddy breakup model of Spalding (1971) is included. This 
model can be written as 

S 0 = C EBU 2 p 0 (2.62) 
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which is based on the idea that the mean reaction rate is 
determined completely by the rate of scale reduction via a 
process of turbulence vortex stretching. Generally, one uses 
the smaller of these two expressions at a given point in the 
flow field. The eddy dissipation rate often slowing chemical 
reaction. This basic model has had several refinements, 
culminating in the Eulerian - Lagrangian theory of turbulent 
combustion (Spalding, 1977) . 

The one-step reaction model of section 2.6.1 (popular 
in the simulation of heat release) may be replaced by a 
slightly more sophisticated kinetics model which allows 
prediction of local mass fractions of hydrocarbon fuel, 
carbon dioxide and nitrogen. Consider the two step 
mechanism: 

C X H Y + <1 + I> <°2 + nN 2> - x CO + l H,0 + (f + J ) nR, 

(2.63) 

X CO + | (°2 + nN 2 ) — X C °2 + 2 nN 

(2.64) 

The mass fractions of all chemical species obey the general 
differential equation. These transport equations are all 
similar and contain terms for convection, diffusion and 
source of a general variable <f> (which contains terms 
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describing the generation and consumption 
they all conform to 


- fv, < r 


3xj 


d<]> 

<t> d*i 


> - S . 


and the equations differ not only in 
diffusive flux but also, and primarily, 
source terms as defined in Table 2.1. 


of <f>) . In fact, 
(2.65) 

their effective 
in their final 
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Table 2.1 Source terms for chemical species 
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where 


r i= <1 + I> w o/ w fu 
r 2 = X w CQ / w fu 

r 3= <V2)W H20 / Wfu 

r 4 - °.5 w 02 /w c0 

r 5 = W C0 2 / W C0 


Further, the diffusion coefficient for all species are 
related to the turbulent viscosity via the turbulent prandtl 
number. For the last entries in Table 2.1, a single values 0 
with a zero source term and values 0 and 1 in the air and 
fuel stream, respectively can provide the solutions for <^ A , 

' i>Q and <f>Q via the following relationships 

0= ^A — I — ^A' a i r = H^R 'a ir = $£ ^P'air (2.67) 

^A'fu “ “^A'air ^B'fu “ ^B'air ^C'fu “ ^C'air 

Using Equation (2.66) and assuming that the mass fraction of 
fuel in the primary stream is 1 and the mass fraction of 
oxidant in the secondary stream is R, we have 

®ox = r 4 m C0 + R(1_0) + ( r i +r 2 r 4 ) (ltl fu" 0) (2.68) 
m C0 2 = (-™£u + 6 > " (2.69) 
m H 2 0 = r 3 ( " m fu + 6) ( 2 ‘ 7 °) 
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This scheme involves the solution of the conservation 
equations for unburnt fuel, CO and mixture fraction. The 
reaction rates of fuel and carbon monoxide are given by the 
following relations 

s fu - - min [ c R1 p m (u * , ^ p‘in f ^m 0 e x exp(-El/Rf) ] (2.71) 

S CQ = - min [ C R2 p m CQ i /»fco m ox”h 2 O exp (-E 2 / R ' i ) 1 (2 ‘ 72) 

where C R1 , C R2 are the eddy breakup constants and k-j^ , k 2 
are the pre-exponential constants for Arrhenius rate 
equations. 

Many chemical reaction phenomena posses stiff kinetics; 
that is, they involve a multitude of species with widely 
different reaction rate coefficient. Their equations have to 
be solved along with the more usual fluid dynamics 
equations. Special computational techniques for handling the 
source terms are required. One often solving the associated 
time - dependent problem even when the steady state solution 
is required. Of course, different time steps are chosen for 
the fluid dynamics and chemical kinetic portions of the 
problems. An efficient technique for steady state solution 
of stiff kinetic problems has been developed by Pratt (1976) 


66 



and shows considerable promise for practical engineering 
situations. 

2.6.3 MODEL FOR PREMIXED COMBUSTION 

In many combustion systems, fuel and oxidant are 
premixed prior to reaction. Premixed flames, spreading 
downstream of flame-holders, have been the subject of many 
previous investigations concerned with the development of 
ram jets and afterburners. In contrast to unpremixed flames 
which can be analysed with fast or finite-rate chemistry, 
premixed turbulent flames are kinetically controlled and the 
rate of flame propagation, called the burning velocity, is 
dependent upon chemical composition and rate of chemical 
reaction. Completely premixed flames are seldom found in 
practice for reasons of safety, (for example flashback and 
blow off) and stability. For this reason, and may be these 
type of flames occur less frequently in practical systems, 
they appear to have received less attention than have 
unpremixed flames. It is only in recent years, however, that 
detailed attempts to understand the flow characteristics, 
with the aid of local measurements, have been reported 
(Stevenson et al. ,1983). A detailed discussion of premixed 
turbulent flames is given by Bray (1980) and a review of 
various methods developed for calculating premixed turbulent 
flames is given by Jones and Whitelaw (1982). 
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In general, most work on premixed flames has assumed 
that combustion can be characterized by a global single step 
reaction of the type shown by equation (2.41). In terms of 
the three - component mixture undergoing a simple one step 
reaction, partial differential equations are set up for m fu , 
mixture fraction and stagnation enthalpy, with the first of 
these for mass fraction of fuel requiring the specification 
of its source, the mean formation rate of fuel. In earlier 
sections the Arrhenius expression was introduced, but 
clearly its mean value should not be calculated merely from 
mean values of components in its expression. If fluctuations 
are neglected and the rate evaluated in terms of the mean 
values of temperature and mass fraction, the result can be 
in error by typically one order of magnitude and will 
exhibit a strong dependence on temperature, pressure and 
mixture strength (Jones and Whitelaw, 1982) . In practice, 
experimental results for premixed turbulent flames are only 
weakly dependent on mass fractions, pressure and 
temperature. This fact led to the turbulence dominated mean 
reaction expression called the eddy breakup reaction model 
described earlier. It is based on the idea that the mean 
reaction rate is determined solely by the rate of scale 
reduction via a process of turbulence vortex stretching. The 
model, thus takes no explicit account of chemical kinetics 
and relates to combustion which is entirely mixing 
controlled (Jones and Whitelaw, 1982) . In this situation, it 
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has been shown to be in good accord with the available 
evidence for premixed flames. The model has also been used 
occasionally for non-premixed flames, but it is generally 
inappropriate in this situation and Pope (1977) has shown 
that it does not necessarily provide unique solution. A more 
logical basis for expressions of the eddy breakup type is 
provided by the Bray-Moss (1977) model for premixed flames. 
The Bray-Moss model for premixed combustion has been 
extended by Libby and Bray (1980a, 1980b), who utilize the 
concept of laminar flamelets in order to derive models for 
turbulent transport and dissipation processes in one 
dimensional planer flames. 

For two-step reaction rate, the solution of the 
conservation equations for unburnt fuel and carbon dioxide 
C0 2 (or carbon monoxide CO) is involved . The mass fraction 
of intermidiate species are determined by the following 
algebraic relations 


m co = Z co' r 2 m fu ' ( m co 2 >/ r 5 

m H 2 o = Z H 2 o" r 3 m fu 

%= Z o 2 + r l m fu “ ( r 4 m co)/ r 5 


(2.73) 

(2.74) 

(2.75) 


where 

Z ccT Vfu.p 


(2.76) 
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°2'P 


r l m fu,p 


(2.77) 
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and m fU/P and ™ o2/ p are mass fractions of fuel and oxidant 
in primary stream, respectively. The mean reaction rates are 
calculated from the minimum of the Arrhenius reaction rates 
and the eddy breakup model similar to (2.71) and (2.72). 



CHAPTER 3 


NUMERICAL PROCEDURE 

The governing partial differential equations described 
in the preceeding chapter are coupled and non-linear. At the 
present time, it is not possible to obtain analytical 
solutions to these equations and numerical technique has to 
be used. A computer program has been developed to solve 
axisymmetric elliptic partial differential equations through 
an iterative procedure based on an integral control volume 
analysis with hybrid upwind finite differencing or quadratic 
upwind differencing and staggered grids. The logic behind 
this program is briefly described below. 


3.1 GRID AND ITERATION SEQUENCE 

Before integrating the standard equations over the flow 
domain, a satisfactory grid is required. For hybrid 
differencing purposes the program uses a staggered grid in 
which the velocities are evaluated at the boundaries of 
scalar variable (P,k,e, etc) cells. Hence separate grids 
define the locations of the U- and V-velocities. A portion 
of these three grids is shown in Figure 3.1 . The solution 
domain is arranged so that the outer surfaces of the 
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boundary scalar cell coincide with the physical boundaries 
of the flow field (Figure 3.2). 

The iteration sequence consists of two major parts: the 
SIMPLE solution algorithm of Patankar and Spalding (1972) , 
and the solution of transport equations for the turbulence 
quantities. The SIMPLE algorithm solves for a fixed pressure 
field, by line iteration, sets of difference equations for 
the x and r momentum equations. After each such sweep over 
the solution domain, adjustments are made to the pressure 
field to satisfy continuity along each line of cells. These 
adjustments in turn destrory the compliance of the velocity 
and pressure field with the momentum equations. Transport 
equations for the turbulence quantities are then solved 
using the calculated velocity field . Finally, closure is 
achieved by evaluating the Reynolds stresses using the new 
mean velocity, turbulent kinetic energy, and dissipation 
rate fields. Iterations are carried out until the momentum 
and continuity equations are simultaneously satisfied to the 
required degree of accuracy. 


3.2 FALSE DIFFUSION 

It is a well-known fact that all upwind scheme, 
although numerically very stable, introduce false diffusion 
into the formulation. It is this false diffusion that acts 
to stablize the numerical solution as it becomes larger than 
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the physical diffusion for P e >2. The exact value of false 
diffusion for one dimensional convection-diffusion equation 
is (Patankar, 1980) : 


I u | ax | P e | r 

r = = (3.i) 

X, f 2 2 


This is not serious when the streamlines are aligned with 
the grid lines, because, in this case, false diffusion 
occures in the direction of the velocity vector and is of 
little consequence for larger P e numbers. However, when the 
streamlines are at an angle to the grid lines, false 
diffusion in the x and r directions combine to give a 
diffusive flux normal to the velocity vector which can lead 
to large errors when the variable considered has a 
significant gradient in this normal direction. That the 
upwind scheme must lead to errors in such a situation 
becomes clear immediately from Figure 3 . 3 which shows that 
it is not the value $ u that is convected across the west 
face but the value $ w _ sg indicated in this figure. An 
approximate expression for the false diffusion coefficient 
for a two dimensional situation has been givened by de Vahl 
Davis and Mallinson (1972). It is 
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p U Ax Ar sin(20) 

4 (Ar sin 3 0+Ax cos 3 6 ) 


(3.2) 


r 


false 


where U is the resultant velocity, and 6 is its angle to the 
x axis. It is seen from this equation that the false 
diffusion is most serous when the flow direction makes an 

o 

angle of 45 with the grid lines. False diffusion can be 
reduced by using smaller grid size and by orienting the grid 
such that the grid lines are more or less aligned with the 
flow direction or using higher order schemes. 

The recently developed skew-upwind and quadratic- 
interpolation schemes have been shown to perform better in 
all situations where streamline skewness causes numerical 
diffusion. Both schemes produce similar results. According 
to Leschziner (1980) , the quadratic scheme requires no more 
computing time than the hybrid scheme, but the skew-upwind 
scheme needs about 50 percent more time. For this reason and 
also because it is easier to implement, the quadratic scheme 
became more popular and is explained in the next section. 

3 . 3 QUADRATIC UPWIND DIFFERENCING SCHEME 

Upwind differencing scheme, although very stable, is 
only first-order accurate and suffers from false diffusion. 
A promising new technique proposed by Leonard (1979), which 
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is known as QUICK (Quadratic Upstream Interpolation for 
Convective Kinematic) , aims at combining the relatively high 
accuracy of the central difference scheme with the stability 
of the upwind scheme. QUICK uses a three-point upstream- 
wighted quadratic interpolation for each cell's wall value, 
two located on either side of the face and the third being 
the next node in the upstream direction. 

Using quadratic interpolation for the dependent 

variable it can be shown that the finite difference form 

of any convective-diffusive equation is reduced to 

A $ =A $ + A$ + A$ + A<& + 

PP EE WU NN SS 

A $ + A $ + A $ + A $ 

EE EE UU UW NN NN SS SS 

( 3 . 3 ) 

The derivation of $ w is given in detail below and the 
derivation of $ at the other faces are similiar. 

A general second degree polynomial for the calculation 
of $ w is 

* w = C 0 + c ! x + c 2 x2 < 3 * 4 ) 

The mathematically equivalent but more convenient quadratic 
form is 
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(3.5) 


* = c o + c i( x-x p) + c 2 (x-x p ) (x-x w ) 

where the constants C Q/ C-^and C 2 are determined from the 
neighbouring points. Consider the situation for positive 
axial velocity in which case calculation of $ w is biased 
toward the upwind node at , then : 

C = $ , (3.6) 

0 P 



(3.7) 


$ WW " $ P % “ $ P 1 

c =( )( ) (3.8) 

2 X WW X P X W - Xp - x w 


and for the case when the axial velocity is negative, the 
coefficient in Eq.(3.5) are calculated using the values of $ 
at the grid node W, P, and E so that: 


c = * / (3.9) 

0 P 



(3.10) 
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1 


) 


(3.11) 


C 

2 



*W ~ *P 
X W X P X E ” X W 


The quadratic interpolation scheme is therefore third-order 
accurate; requires 9-point and must take the sign of all 
face velocities into account. This algorithm was extended to 
two dimension by applying a similiar procedure in the cross- 
stream direction. Like the central-difference scheme, the 
quadratic scheme does not produce diffusion-type truncation 
error, but it may suffer from some unboundedness. Han et al. 
(1981) further report that, unlike the upwind scheme, the 
quadratic scheme is not unconditionally stable and may 
require certain special measures to produce stability. 
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Solution domain used in backward- facing 
step flow computations 
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Figure 3.3 Quadratic upstream-weighted difference scheme 
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CHAPTER 4 


MODEL EVALUATIONS 

This chapter presents the results of the comparisons of 
the various models for isothermal, constant and variable 
density turbulent recirculating flows. Predictions obtained 
with each of the turbulence models are discussed and, when 
available, compared with experimental data. The goals of 
this study require a careful selection of the test cases. 
They have to provide reliable mean flow and turbulence data 
in the recirculation and recovery regions. A detailed 
specification of the flow parameters in the upstream region 
is also essential since these are used as inlet conditions 
start the computations. However, the most common 

shortcoming is found to be the lack of well-defined inlet 
boundary conditions. The importance of these quantities has 
not been fully appreciated (Sturgess et al . ,1983). 

Incorrect specification of inlet turbulence quantities could 
have a further adverse effect if the flow is a reacting one. 
This is because the eddy break-up combustion model relates 
the reaction rate of the fuel to the eddy lifetime, k/e. An 
incorrect specification • of k and 6 at the inlet therefore 
results in an incorrect density field. Since density is 
strongly coupled to the mean flow field, an incorrect 
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density changes the mean flow field completely. Finally, all 
data sets have to meet the criteria set by Eaton and 
Johnson (1980) for acceptable measurements in terms of 
adequate experimental facility, appropriate instrumentation, 
and agreement with generally accepted flow trend. 

After reviewing some of the available data sets the 
following three are chosen as the test cases: 

1. Johnson and Bennett (1981); confined coaxial 
non-swirling suddenly-expanded jets 

2. Roback and Johnson (1983); confined coaxial 
swirling suddenly-expanded jets 

3. Brum and Samuelsen (1982); C0 2 into air 

(three coaxial streams with middle stream swirling) 

Only a portion of the predictions and comparisons with 
measurements are presented. These are chosen to illustrate 
specific points concerning the turbulence models and their 
performance. 
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4.1 THE BASIC EXPERIMENTAL FLOW FIELDS 


The cases of Johnson and Bennett (1981) , and Roback and 
Johnson (1983) are selected for comparison with the 
calculated mean and turbulence quantities in a constant 
density flow. The experiments are carried out in a water rig 
that resembles closely a gas turbine engine (can) combustor 
(Figure 4.1). The flow conditions selected for this study 
have Reynolds numbers of 15900 and 47500 for the inner and 
annular streams, respectively. These Reynolds numbers are 
factors of 5 to 20 greater than the transitional Reynolds 
number range occurring in aircraft gas turbines. Mean 
velocity and turbulence fluctuations were measured with a 
Laser Doppler Velocimeter ( LDV) . The scalar transport 
measurements were made using a Laser-Induced Fluorescence 
(LIF) technique. Fluorescein dye was introduced as a tracer 
into the water flow of the central jet to give constant 
density, variable species, mixing of the two concentric 
jets. This permitted the modelling of turbulent mass 
transport to be investigated. Sketches of the flow regions 
occurring for the swirling and non-swirling flow conditions 
are shown in Figure 4.2. For the non-swirling flow, four 
major shear flow regions can be identified. These are: a 

wake region immediately downstream of the inner jet inlet 
duct, a shear region further downstream between the inner 
and annular jets, a recirculation zone, and a reattachment 
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region. A three-dimensional counter-rotating corner eddy is 
also observed and extends only a short distance, less than 
one step height downstream of the step. The large primary 
eddy occupies the rest of the recirculation zone. For 
swirling flow, a second recirculation cell occurs along the 
centerline and the length of the annular recirculation cell 
is decreased. One other difference between the swirling and 
non-swirling experiments should be noted. This is the 
observation that, once reattachment of the annular jet 
occurred, the axial velocity profile tended to flatten with 
momentum transport from the outside inward for the swirling 
flow case whereas the velocity profile tended to flatten 
with momentum transport from the inside outward for the 
nonswirling flow case. 

For variable-density flow analysis, the isothermal, 
non-reacting flow inside a dilute swirl combustor provided 
by Brum and Samuelsen (1982) is selected. The dilute swirl 
combustor is designed with a cone— annular gas injector, a 
swirling stream and a non— swirling stream of dilute air. The 
complex model laboratory flow combustor (Figure 4.3) has an 
aerodynamically controlled, swirl— stabilized recirculation 
zone. It consists of an 80 mm I.D. cylindrical stainless 
steel tube 50 cm long with rectangular optical windows 
mounted vertically on either side of the combustor tube. 
These flat windows provide for clear optical access 
necessary for laser anemometry measurements. A set of swirl 
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vanes (57 min O.D.) is concentrically located within the 
combustor tube around a 19 mm O.D. centrally positioned 
fuel/C0 2 delivery tube. The combustor was operated at 
atmospheric pressure with an overall equivalence ratio of 
0.2 and a bulk reference velocity of 15 m/s. 

4.2 EVALUATION OF DIFFERENT CLOSURE MODELS 
FOR COMBUSTOR FLOW CALCULATIONS WITHOUT SWIRL 

In this section the comparison is made of the 
performance of the k-e closure, algebraic stress closures 
and Reynolds stress closure in calculating the combustor 
flows. Effects of different pressure-strain models and 
turbulent diffusion models are also analysed. Predictions of 
the mean and turbulence quantities using high and low 
Reynolds number models are compared. Finally, the efficient 
models for turbulent momentum exchange in combustor flow 
calculations are identified. 

4.2.1 THE k-e MODEL RESULTS 

The mean and turbulence quantities predicted by the 
standard k-e model and the measured values are illustrated 
in Figures 4.4 - 4.6. The model constants used in this 

calculation are listed in Table 4.1. The overall agreement 
is good and the calculated velocity and maxima and their 
locations closely reproduce the experimental results. The 
recirculation zone and the reattachment length has been 
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reasonably predicted (Tables 4. 2-4. 3). There are however, 
some discrepancies in the region of the center line with the 
velocities being under-predicted by about 5% - 30% at the 
downstream locations. The change in the axial velocity 
profile from x = 13 mm to successive locations downstream 
document the development of the various shear regions within 
the combustor. The shear layer between jets occurred between 
x=51 mm and 203 mm. In this shear region, the annular jet 
flow accelerated the inner jet flow but the predicted result 
does not show such a behavior and the centerline velocity is 
contineously decreased (Figure 4.7). The under-prediction 
for a centerline velocity stems from the incorrect 
representation of the turbulent diffusion process. As shown 
by Ribeiro (1976) , the radial normal stress is particularly 
important in the upstream region and, as a consequence, the 
isotropic viscosity hypothesis is inadequate. 

Comparison of the predicted turbulent kinetic energy 
with data are shown in Figure 4.5. The predicted k values 
are in good agreement with the data near the inlet plane. 
However, the k-e model under-predicts the k values further 
downstream near the centerline. The discrepancy near the 
centerline is due to: (1) under-prediction of the 
corresponding maximum in the mean velocity, (2) inadequacy 
of the eddy viscosity hypothesis, and (3) over-prediction of 
turbulent kinetic energy dissipation rate, e. The inlet 
distribution of length scale has a significant effect on the 
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e distribution in the upstream region where the diffusion 
terms have a large effect on the mean velocity; it has no 
significant effect in the downstream region where the 
diffusion effects are smaller. The k-e model tends to 
connect the dissipation rate too strongly to the local mean 
velocity field (Habib and Whitelaw, 1979) which is 

inappropriate since the dissipation occurs in the finest 
scales of motion and these do not reflect the local mean 
strain field. These connection tends to increase the local 
level of e. The predicted turbulent shear stress ,uv, 
profiles are also in very good agreement with the data 
(Figure 4.6) . At x = 13 mm, the negative peak value for uv 
near the axis corresponds to the shear layer between the 
inner and annular streams and the positive peak corresponds 
to the shear layer associated with the pipe expansion. 


4.2.2 THE ALGEBRAIC STRESS MODEL RESULTS 

Prediction results for non-swirling flow using n^j 
model given by eq. (2.13) are shown in Figures 4.8-4.13, 
while the comparisons of different n^j models given by 
eqs. (2 . 13 , 2 . 14) , (2.13,2.15) and (2.13,2.16) are presented 

in Figures 4.14-15. The various models used in this 
investigation and the recommended values for the constants 
appearing in these models are presented in Tables 4.4-4. 6. 
The details of these calculations have been discussed by 
Nikjooy et al. (1985). 
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For a fixed pressure-strain model, the calculated mean 
and turbulence fields obtained from the three different ASM 
are essentially identical and are in good agreement with 
measurements. Some slight differences can be identified; 
however, they are not substantial enough to lead one to 
conclude that one ASM is better or worse than another. In 
general, the flow is very well predicted in the near and far 
field but not so well calculated in the mid-field. This is 
especially true of the flow near the combustor core, where 
the calculated mean axial velocity and turbulent normal 
stresses are consistently lowered than the measurements 
(Figures 4.8,4.11-4.13). The recirculation zone is very well 
predicted by all ASM. One surprising result of the present 
study is given by the calculations of the k-e closure 
(Figures 4.8-4.10). Essentially the same calculated 
distributions of U, k and uv are obtained when the k-e 
closure is used instead of the ASM. Indeed, the k-e closure 
gives a k distribution that is in as good an agreement with 
measurement as that given by ASM/E (Figure 4.9). 

A second interesting result of the present study can be 
found in the comparison of the different n^j models. Here, 
the effects of the four n^j models are discussed. Instead of 
giving rise to improved correlations with measurements as 
anticipated, some of the calculations are actually in worse 
agreement with measured data for ASM/MY and ASM/R (Figures 
4.14,4.15). As for the ASM/E, very little differences are 



noted for the four n^j considered, except in the prediction 
of the reattachment length. Also, the choice of the mean 
strain model has a great effect on the calculated turbulence 
field. The results shown in Figures 4.14 and 4.15 for the 
11^ given by (2.13,2.14), (2.13,2.15) and (2.13,2.16), 

respectively, are selected to illustrate the points made 
above. In general, the comparisons for other stream 
locations are quite similar to those shown in Figures 4.14 
and 4.15. Depending on the non-equilibrium ASM used for 
closure of the flow equations, n^j, 2 modelling given by 

(2.14) or (2.15) will lead to prediction of double peaks 
behavior for the turbulence field (Figure 4.14), which are 
not observed in the experimental flow. Such behavior is not 
found in ASM/E calculations, though. One reason for this 
could be due to the fact that the constant C 2 used is not 
suitable for complex turbulent flows. Since C 2 is determined 
for simple turbulent flows in local equilibrium, it would be 
more appropriate for ASM/E than for ASM/MY and ASM/R. 
Another reason could be the incorrect modelling of the mean 
strain part of the pressure-strain terms by (2.14) and 

(2.15) . The results shown in part (c) of Figures 4.10 and 
4.11 illustrate the importance of accounting for the anti- 
symmetric contributions of the mean strain tensor as given 
by (2.16). With this improvement, the anomalies seen in 
parts (a) and (b) of Figures 4.14 and 4.15 disappear with 
the exception of the calculated -Jv* using ASM/R closure. 
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Therefore, it is very important to model the mean strain 
tensor correctly. Either the mean strain effects are not 
modelled at all, as in (2.13), or they should be accounted 
for properly, as in (2.16). 

Another measure of the performance of the different 
turbulence models is in the calculation of the reattachment 
length. This length is determined by locating the point of 
zero shear at the combustor wall. Factors which influence 
the reattachment length may be classified in one of two 
groups: (1) system geometry, (2) upstream conditions. System 

geometry has a significant effect on reattachment length. 
The two primary geometric parameters influencing the sudden 
expansion flow are: (a) d 2 /d 1 , the expansion ratio, where d-^ 

is the diameter at the inlet to the sudden expansion and d 2 
is the diameter of the downstream tube, (b) d 2 /(d 2 -d 1 ) ,the 
aspect ratio. Upstream conditions which have an effect on 
reattachment length in the sudden expansion flow are: (a) 

the inlet flow Reynolds number, (b) centerline turbulence 
level, (c) inlet Mach number. So (1986) has examined the 
effects of these parameters on reattachment length in an 
axisymmetric sudden-expansion flow. Based on his analysis, 
it was reported that the single most important parameter 
that influences the reattachment length is the centerline 
turbulence level. The results of the fifteen different 
calculations are listed in Table 4.3 together with the 
measured reattachment length determined from Johnson and 
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Bennett (1981) . It can be seen that the best estimate of the 
reattachment length is given by the k-e and ASM/E closures, 
while the worst is given by the ASM/R with n^j deduced from 
(2 . 15) . 

An Attempt is made to assess the lack of differences 
shown by the calculations of the three ASM closures (Figures 
4.8-4.13). The calculated turbulent viscosities of the three 
ASM closures are compared with the numerical viscosity 
(Figure 4.16). It can be seen that the turbulent viscosities 
are quite a bit smaller than the numerical viscosity in the 
inlet region of the combustor and are of the same order at 
or near the dividing streamline of the recirculation zone. 
Consequently, calculations in the near field and the 
recirculation region are being masked by the numerical 
viscosity and the errors would propagate through the whole 
field. Unless the numerical viscosity is significantly 
reduced, the improvements afforded by more sophisticate 


turbulence 

closures will not 

be 

realized 

for 

complex 

turbulent 

flow calculations. 

In 

addition 

to 

hybrid 


differencing scheme, the Quadratic Upwind Differencing (QUD) 
scheme of Leonard (1979) was also explored. The results are 
very much similar to the hybrid scheme predictions. It 
should be pointed out that QUD reduces false diffusion, but 
does not eliminate it entirely (Figure 4.17) 
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4.2.3 REYNOLDS STRESS MODEL RESULTS 


Having demonstrated the effect of the pressure-strain 
correlation and compared the k-e model with ASM, attention 
is now turned to the performance of the Reynolds-stress 
model (RSM) . It seems somewhat paradoxical that the more 
advanced types of turbulence model have been least 
successful in the complex flows associated with flow 
recirculation for which (in theory) they have the most to 
offer compared with eddy viscosity based closures. The 
difficulty in clearly demonstrating calculations free of 
numerical errors has hampered the real testing of Reynolds 
stress closures in such flows due to the practice of using 
upwind differencing for discreption of the convection terms. 
This practice is resulted in the phenomenon of numerical 
diffusion and to produce a scheme which is slow to respond 
to grid refinement. Although local selective mesh refinement 
can be used to obtain numerical error-free predictions when 
using eddy viscosity based models, this become totally 
impracticable when considering Reynolds stress transport 
closures where the lack of an eddy viscosity increases the 
local cell Peclet numbers and the associated numerical 
diffusion coef f iecients by orders of magnitude. Unless 
numerically accurate solutions are developed, the 
improvements afforded by advanced turbulence models will not 
be realized for complex recirculating flows. 
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Effects of two different turbulent diffusion models 
(2.18,2.19) on the Reynolds stress closure are investigated 
(Figures 4.18-4.21). The predicted results of the mean 
velocity and turbulent shear stress, uv fields obtained from 
these two models are very similar. However, some differences 
are observed in the prediction of the peak values of the 
normal stresses. In the fully developed region, the model 
developed by Daly and Harlow (2.18) performs better, 
especially near the centerline. Despite the superiority of 
Hanjalic and Launder' s model (2.19) in a theoretical sense, 
it does not appear to give better predictions. 

For a fixed diffusion model (2.19), two different 
pressure-strain models are employed to predict the mean and 
turbulence quantities. In one case, mean strain rate effects 
are not considered and in another, these effects are 
accounted for using the Launder et al. model (2.16). The 
predicted results of mean velocity and turbulent stresses 
are illustrated in Figures 4.22-4.25. It can be seen that, 
both models performed similarly. However, Rotta's model 
(2.13), gives a better prediction of the normal stresses. 

In comparison with other simple closures (k-e and ASM) , 
it can be seen that the major discrepancy arises in the 
prediction of reattachment length (Table 4.3) and centerline 
values (Figures 4.26-4.30). In spite of its higher class of 
sophistication, RSM under-predicts the size of the 
recirculation region. A question arises here on the mass 
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flow rate in the combustor displayed by the measured data; 
this seems larger than in the calculations, which certainly 
conserve mass between inlet and exit from the solution 
domain. This feature is also to be seen in the comparisons 
shown by Mansour et al. (1983) and Hackman et al. (1984), 
where the area under the predicted velocity curves is 
clearly less than under the measured data. It is possible 
that transverse side-wall boundary layer growth in the 
experiment caused acceleration on the centerline of the 
combustor, but it is difficult to do more than speculate. 

It is speculated by McGuirk et al. (1985) that the k-e 
deficiency is caused by the model's inability to represent 
the normal stress-normal strain production terms properly. 
This feature still remains in the RSM calculations, as is 
shown by the distribution of the two contributions to the 
production of k from shear and normal stress terms in the 
profiles drawn in Figure 4.31 at an axial station near 
reattachment. Although positive values of the normal stress 
term are obtained near the axis, they represent only about 
10% of the maximum shear stress production. Estimate made 
from the measurements of recirculating flow reported by 
Taylor and Whitelaw (1984) show that the maximum values 
should be of the same order, so that even in the RSM 
calculations the normal stress production is too small by an 
order of magnitude. 
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So far computations have been performed by using a high 
Reynolds number model and 51x41 non-uniform grid system. The 
conventional wall law has been adopted to connect adjacent 
grid lines to the wall boundary. However, results do not 
show any trace of the secondary circulation observed 
experimentally. The switchover to a fine-grid (51x50) , low- 
Reynolds-number analysis has had a significant effect on the 
predicted flow pattern in this region. An appreciable 
secondary eddy is now formed which extends almost one— step 
height downstream and one-tenth of a step height normal to 
the wall (Figure 4.32). It is because of the thinness of 
this secondary eddy that the wall-function approach failed 
to predict its existence. Comparison of the calculations 
with a high and a low Reynolds number model results show 
that outside the viscosity dominated region, the two models 
performed with just about equal success (Figures 4.33-4.34). 
However, transpiration through the wall and steep 
temperature gradients due to large imposed wall heat fluxes 
or frictional heating are just a few of the cases that may 
cause the near wall region to have more influence on the 
core region. Despite the limitations of the turbulence 
model, the replacement of wall function by a fine-grid 
treatment has allowed a more realistic modelling of the flow 
just downstream of the step. 
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4.2.4 CONCLUSIONS 


Based upon the preceding discussion, Several main 
conclusions can be summarized as follows. 

1. As far as the mean field prediction is concerned, k-c 
closure performs just as well as any ASM closure in the 
calculation of combustor flow. 

2. The reattachment length is reasonably predicted by the k- 
£ closure. 

3. The effects of pressure-strain modelling on the 
calculations are large compared to the effects of different 
ASM closures. 

4. If the mean strain effects on the pressure-strain terms 
are to be accounted for, they should be modeled by the 
Launder et al. (1975) model. 

5. The simplest pressure-strain model (2.13) gives as good a 
correlation with measurements as the model given by (2.13, 
2.16) 

6. The best correlation with measurements for combustor flow 
calculations is given by the ASM/E closure with 11^ 
determined from (2.13). 

7. The best prediction of normal stresses is given by the 
ASM/E closure with n^j determined from (2.13). 

8. The best prediction of shear stress is given by the ASM/E 
closure with n^-s determined from (2.13) or the k-£ closure. 
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9. Turbulent diffusion model developed by Daly and Harlow 
(1970) results in better predictions than the Hanjalic and 
Launder (1972) model. 

10. Reynolds stress models do not seem to yield better 
predictions than the algebraic stress models. 

11. Low-Reynolds-number model prediction of the near wall 
region is more realistic and show the corner recirculation 
region. 


4.3 MODEL EVALUATIONS FOR COMBUSTOR FLOW 
CALCULATIONS WITH SWIRL 

In this section, performance of the k-e closure, 
algebraic stress model (ASM) and Reynolds stress model (RSM) 
are evaluated and compared with data. In the previous 
section, it was concluded that the turbulent diffusion model 
developed by Daly and Harlow (1970) results in better 
predictions than the Hanjalic and Launder (1972) . Therefore, 
this model is employed for RSM calculations. Based on some 
preliminary calculations, it was found that the mean strain 
rate effects have to be accounted for. Otherwise, predicted 
results of the normal stresses show negative values in some 
regions of the flow fields. Therefore, these effects are 
approximated by the Launder et al. model (2.16) 

The data of Roback and Johnson (1983) represent a 
carefully measured turbulent water case with a 30 degree, 
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free swirler. Inasmuch as inlet profiles were not provided, 
Calculations were started at the first measurement location, 
which was just downstream of the expansion plane (5 mm) . 
This approach is unsatisfactory because the input length 
scales of turbulence used to calculate dissipation rate are 
unknown and the effect of area expansion on the flow has not 
been properly accounted for. The ability to predict the 
existence and location of the central recirculation zone in 
this combustor is a good indication of the suitability of 
the model . 

4.3.1 THE k-e MODEL RESULTS 

The centerline plot of the axial velocity is provided 
in Figure 4.35. The k-e model mimics the data trend 
reasonably well along the symmetry axis, however, the rate 
of recovery is under-predicted. The reason is not easily 
resolved. Comparison of the calculated mean axial and 
tangential (azimuthal) velocity profiles with the 
experimental data are presented in Figures 4.36 and 4.37. In 
this case, the k-e model fails to display the size and 
strength of the experimental recirculation zone. The 
predicted axial profile is skewed tward the wall in 
accordance with the data trend, although the maximum value 
decreases slightly with a corresponding increase in the 
near-wall region. At the farthest downstream locations, the 
k-e model produces a decreased axial velocity near the 
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centerline, indicative of its lagging characteristics and 
apparent slow recovery. The k-e model prediction of 
tangential velocity is excellent up to a station of 152 mm 
from the expansion. At subsequent downstream locations, 
calculation of the azimuthal velocity decays prematurely to 
a forced vortex profile, while the experimental data still 
shows a combined vortex profile. Over most of the radius, 
the measured tangential velocity is approximately constant 
and only rapidly approaches zero at the centerline from a 
radius of about 2 cm (r/R= 0.33). Further downstream from 
305 mm (x/R=5) , the calculated flow approximates a forced 
vortex while the measured flow approximates a free vortex in 
the outer region, with a forced vortex core. 

4.3.2 THE ALGEBRAIC STRESS MODEL RESULTS 

For a pressure correlation model given by Launder et 
al.(2.16), three different ASM's are compared. Results show 
that the calculated mean and turbulence fields are very 
similar for all the ASM's and are in good agreement with 
measurements at most of the locations (Figures 4.38-4.41). 
The tangential velocities obtained from the three ASM's are 
also very similar and in very good agreement with the data. 
The algebraic stress models provide profile shapes which are 
remarkably similar to the exhibited data trend, although the 
swirl velocity predictions are reduced in magnitude from the 
measurements. The discrepancy in tangential velocity 
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profiles is probably generated in part by incorrect inlet 
radial velocity profile. The measurements of mean axial and 
radial velocities at downstream stations indicate that the 
tangential velocity profile at 405 mm is in equilibrium. A 
free vortex is not a stable profile, and viscosity effects 
at the core modify the profile to the stable forced vortex 
form locally, as indicated at 405 mm. Therefore, the 
discrepancy in tangential profiles could also be due in part 
to an inadequate calculation of eddy viscosity across the 
radius . 

Figure 4.40 displays the fluctuating axial velocity at 
various locations. All the algebraic stress predictions are 
similar to the exhibited data trend with a slight under 
prediction of stress magnitudes. This may be due to the 
modelling of the terras in the pressure-strain correlation 
and the constants used. Consequently, it cannot be implied 
that the algebraic stress model exhibits a superior 
prediction capability. The lack of agreement between the 
predictions and the data may be attributed to experimental 
error, boundary conditions, numerical diffusion, and 
oscillatory phenomena. However, it is difficult to separate 
the relative effects and apportion to each its respective 
contribution to error. 
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4.3.3 REYNOLDS STRESS MODEL RESULTS 


Results of the mean and turbulence field obtained from 
the k-e model and the ASM were compared with the Reynolds 
stress model (Figures 4.42-4.45). A wide disparity exists 
between the models prediction of the axial velocity near the 
centerline. It is evident that there are substantial 
differences in the capability of the various models to 
promote or hinder formation of the recirculation zone 
relative to the k-e model . The comparison of the calculated 
and measured locations of the forward and rear stagnation 
points along the centerline for five different models is 
presented in Table 4.7. It is clear that the proper 
turbulence model is dependent on the location within the 
flow field. With regard to the comparison between 
measurements and calculations, the predictions by RSM seem 
to be slightly better than ASM for tangential velocity, 
while those by k-e model are closer to experimental data for 
the central recirculation. An examination of the calculated 
Reynolds stresses indicates that the relative performance of 
the model is strongly dependent on the flow region. 

4.3.4 CONCLUSIONS 

Based upon the preceding discussion, The main 
conclusions can be summarized as follows. 
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1. As far as the mean field prediction is concerned, 
algebraic stress model (ASM) and Reynolds stress model (RSM) 
perform better than the k-c closure in predicting the 
tangential velocity profiles. 

2. The mean strain effects on the pressure-strain terms must 
be accounted for. 

3. The central recirculation region is over-predicted by the 
ASM and the RSM; however, it is under-predicted by the k-e 
closure. 

4. The Reynolds stress model does not yield any better 
predictions compared to the algebraic stress model 
calculations . 

5. Inlet conditions are the most important factor in 
determining the location, size, and the strength of the 
central recirculation region. 


4.4 SCALAR TRANSPORT MODELLING AND COMPARISON 

In the present study, the experiments of Johnson and 
Bennett (1981) and Roback and Johnson (1983) are chosen to 
compare with the calculation results. These experiments are 
selected because of their unique turbulent mass transport 
measurements . 
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4.4.1 NON-SWIRLING FLOW CALCULATIONS 


In the present study, the mean velocities and Reynolds 
stresses predicted by Reynolds stress closure are used to 
calculate the mean mixture fraction, 0, the RMS of the 
fluctuating mixture fraction and turbulent mass flux for 
non-swirling flow (Figures 4.46-4.50). The various models 
used for turbulent scalar flux and the recommended values 
for the constants appearing in these models are listed in 
Table 4.8-4.10. In the present study, firstly, the effect of 
the two algebraic flux closures (2.36,2.37) with two 
different pressure-scalar gradient models on combustor flow 
mass transfer are analysed (Figures 4.46-4.48). Secondly, 
the components of the scalar flux are obtained directly from 
solution of their respective modelled transport equations. 
The effect of pressure-scalar gradient models are also 
investigated. All the model predictions are essentially 
identical and are in good agreement with measurements. Some 
slight differences can be identified, however, they are not 
substaintial . In the case of mean concentrations, the 
calculated profiles at x=51 mm, x=102 mm and x=3 05 mm, are 
in close agreement with those measured. Between x=152 mm and 
x=254 mm where the transport in the axial direction is 
comparable to or larger than the radial, some discrepancies 
near the axis of the tube are seen. Here, the turbulent 
diffusion rates are anisotropic. Although this effect has 
been considered in the modelling, the model did not respond 
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effectively in the non-isotropic region of the flow field. 
The shapes and thickness of the concentration profiles are 
in reasonable agreement with measured values, however, the 
calculated concentrations at and in the vicinity of the 
center line are over-predicted by around 2 0% . This is 

consistent and at least in part a consequence of the low 
mean velocities calculated in this region. The error is 
probably due to experimental errors, because the predicted 
levels of the mixture fraction are required to conserve the 
mass flux of dye. The fluctuation profiles show that the 
models under-predict the peak values almost at all the 
locations. These predictions suggest that the constant 
associated with the dissipation of fluctuations should be 
lower than the present choice in some portions of the flow 
field. 

The results of the transport equations for turbulent 
scalar fluxes (Figures 4.49-4.50) show that the predicted 
mean concentration are in better agreement with the data 
although the turbulent mass flux has been over-predicted 
near the centerline. It seems reasonable to increase the 
diffusion rate of the turbulent mass flux through its 
constant. 

Finally, it should be pointed out that Mellor and 
Yamada's model gives a linear set of algebraic equations for 
uTUj. and 5u. . This is much more easy to solve and less time 
consuming than the set of non-linear algebraic equations 
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obtained from Rodi's model (2.33) ,and Launder and Gibson's 
model (2.37) . 

4.4.2 SWIRLING FLOW CALCULATIONS 

For swirling flows, effects of two pressure-scalar 
gradient correlation models (2.22,2.23) on two Algebraic 
Flux Models (AFM) are also investigated (Figures 4.51-4.53). 
The calculated mean and turbulence fields obtained from 
these closure models are very similar. Some minor 
differences are observed. However, they are not significant. 
The agreement between the predicted and the measured mean 
concentration profiles at most of the locations are 
excellent, except near the axis at the central recirculation 
zone (x=51 mm) . The discrepancy could be due in part to an 
inadequate calculation of the centerline velocity. The 
radial turbulent mass transport rate profiles indicate that 
radial transport occurred at axial locations of x= 13, 25, 

50 mm and essentially zero mass turbulent transport between 
x=102 mm to x=203 mm. The peak rates for x=13 mm and 25 mm 
correspond to the location of the interface between the 
inner jet and annular jet, and the location where the axial 
velocities were negative, respectively. At these two 
stations the turbulent mass rates have been under-predicted 
which could be the results of the inlet condition or 
velocity field. The discrepancy between the model 
predictions and the measurement is more clear at 51 mm. 
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There is considerable scatter in the measurements, however, 
this location is near the upstream end of the central 
recirculation region where the flow may not be axisymmetric. 
Finally, the ^-probability density function for 
concentration at location x = 2 5 mm is calculated from the 
known mean and variance. The results obtained are very 
similar and in close agreement with the measured profiles 
(Figure 4.54) . 

4.4.3 CONCLUSIONS 

The main conclusions emerged from the preceding 
discussion can be summarized as follows: 

1. The two algebraic flux model predictions are very 
similar. 

2. The effect of mean strain modelling on the pressure- 
scalar gradient terms is not significant. 

3 . The transport model for scalar fluxes does not have any 
advantage over the algebraic flux model. 

4. The /^-distribution is sufficient to give the correct 
behaviour of the instantaneous scalar distribution. 
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4.5 EXTENSION OF CONSTANT DENSITY MODELLING TO 


VARIABLE DENSITY CALCULATIONS 

The preceding models are extended to the variable- 
density swirling flow conditions studied by Brum and 
Samuelsen (1982) . For this case, mean velocity and 
turbulence intensity measurements were provided in the flow 
field for the axial and tangential velocity components. 
However, inlet velocity profiles were not measured. 
Consequently, the axial velocity profile used to initiate 
computations was regarded to an idealized, turbulent, 
annular pipe flow distribution. The tangential mean velocity 
component and turbulence quantities were taken from the 
downstream measured profiles. The turbulence length scale at 
inlet was assumed to be 1.6 mm. It may be argued that the 
independent specification of the velocity distributions for 
separate jet streams, the use of idealized profiles, or the 
use of profile shapes from the measurements is not entirely 
justifiable. Nevertheless, the predictions are still useful 
in demonstrating some of the differences between the various 
turbulence models. The injection velocity of propane/C0 2 was 
estimated based upon the mass flow and nozzle area. It is 
quite high (28.1 m/s) and such velocities would tend to 
eliminate or change the position of the recirculation 
region. Therefore, accurate inlet injection velocity 
measurements are required to correctly predict the flow 
field in the recirculation zone. 
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4.5.1 k-e MODEL VERSUS ALGEBRAIC STRESS MODEL 

Computations are started from the swirler exit, 
extending 300 mm downstream. A non-uniform grid of 51x50 
with more nodes concentrated from the centerline to the 
nozzle opening is used. Axial grid spacing starts with 0.5 
mm at the inlet and expands geometrically to the exit. A 
total of four different sets of calculations are performed 
for the non-reacting case. The calculations and their 
comparisons with measurement for n^j model given by 
(2.13,16) are presented in Figures 4.55-4.58. The results of 
the three ASM predictions of mean and turbulence fields are 
very similar. In the outer flow region, the predicted axial 
velocities are in agreement with the data; however, the 
calculations reveal substantial differences between the 
ASM's and the k-e model near the centerline. The differences 
appear to be more significant for the present case than the 
previous ones, because of the large impact on the primary 
jet by the surrounding, swirling flow field. All the models 
induce a recirculation zone along the symmetry axis but 
maximum flow reversal velocities are smaller than the 
measurement. A comparison between predicted and measured 
swirl velocity profiles show that the algebraic stress model 
promotes fairly rapid decay of the mean velocities relative 
to the k-e model. The use of the algebraic stress model, 
which are destabilized with forced vortex flow, may not be 
sound. The predicted turbulence profiles demonstrate a 
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pronounced relaminarization effect, reducing the turbulence 
intensities to unrealistically low values. The discrepancy 
between the prediction and the measurement is probably 
generated in part by incorrect inlet condition and to some 
extent by using constant-density models for variable-density 
flow. However, this is minor, since the variations in 
density flow field are not significant enough to affect the 
momentum field (air flow rate/C0 2 flow rate = 117.6). It is 
recognized that by selecting only three examples, it is 
difficult to generalize about the performance of these 
models under all conditions. It should be noted that this 
particular set of constants for the ASM cannot be lightly 
dismissed. The unusual characteristics adopted by the 
predictions prompts the inquiry as to the necessity of the 
full convection and diffusion terms for strongly swirling 
flow. 

4.5.2 CONCLUSIONS 


The previously discussed comparisons have demonstrated 
the relative merits of various turbulence models for the 
calculation of swirling, recirculating flows. It should be 
noted that the effectiveness of turbulence model predictions 
may be obscured to some extent by competing factors such as: 
inlet and boundary conditions, oscillatory phenomena, and 
numerical scheme. A significant contribution from any of the 
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aforementioned factors tends to invalidate conclusions 
regarding the superiority or inferiority of a given 
turbulence model . Even the most advanced turbulence model 
cannot compensate for inadequacy in this area. 

Although little confidence is expressed in the ability 
of current generation turbulence models to simulate swirling 
flow aerodynamics, the k-e model performs competitively. 
None of the ASM's could satisfactorily predict the Reynolds 
stresses, but alteration of the constants in pressure-strain 
model is a viable option to improve the capability of the 
ASM's for strongly swirling flows. The derivation and 
validation of higher-order closure schemes hold the greatest 
potential for turbulence model improvement for strongly 
swirling flows. In view of the current numerics, the 
additional computational time associated with higher order 
closure, and the lack of a well validated turbulence model, 
the k-e model remains the model of choice and should be used 
to calculate reactive flows in practical combustors. 
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Table 4 . 1 Values of constants in k-e model . 


c 

C .1 

C s 2 

° k 

a 

e 

0.09 

1.44 

1.92 

H 

• 

O 

1.3 


Ill 


Table 4.2 Various models used in Reynolds stress 
and algebraic stress closures 


MODEL 

MODEL NUMBER 

EQUATION 

NUMBER 

REDISTRIBUTION 

1 

(2 . 13) 


2 

(2.13) + (2.14) 


3 

(2.13) + (2.15) 


4 

(2.13) + (2.16) 

DIFFUSION 

5 

(2.18) 


6 

(2.19) 

DISSIPATION 

7 

(2.24) 


8 

(2.40) 
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Table 4 . 3 Comparison of calculated and measured 
reattachment lengths in mm. 


MODEL NUMBER 




1 

2 

3 

4 

k-e 

ASM/E 

258 

245 

242 

227 

246 

ASM/MY 


245 

254 

222 

246 

ASM/R 


240 

229 

158 

201 


RSM 

DATA 254 


113 



Table 4.4 Various algebraic stress models (ASM) 
used in this investigation 


MODEL 

MODEL NUMBER 

EQUATION 

NUMBER 

ASM/E 

9 

(2.28) + (model-1) 


10 

(2.28) + (model-2) 


11 

(2.28) + (model-3) 


12 

(2.28) + (model-4) 

ASM/MY 

13 

(2.31) + (model-1) 


14 

(2.31) + (model-2) 


15 

(2.31) + (model-3) 


16 

(2.31) + (model-4) 

ASM/R 

17 

(2.34) + (model-1) 


18 

(2.34) + (model-2) 


19 

(2.34) + (model-3) 


20 

(2.34) + (model-4) 
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Table 4.5 Various Reynolds stress models (RSM) 
used in this investigation 


MODEL 

MODEL NUMBER 

EQUATION 

NUMBER 

high-Reynolds 
number model 

21 

(2.7) + (model-1) + 
(model-5) + (model-7) 


22 

(2.7) + (model-4) + 
(model-5) + (model-7) 


23 

(2.7) + (model-4) + 
(model-6) + (model-7) 

low-Reynolds 
number model 

24 

(E.1-E.6) + (model-4) 
(model-5) + (model-8) 


25 

(E.1-E.6) + (model-4) 
(model-6) + (model-8) 
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Table 4 . 6 Values of constants used in Reynolds stress 


and 

algebraic stress 

closures 


MODEL 

CONSTANT 

ASSIGNED 

VALUES 

MODEL 

NUMBER 

REDISTRIBUTION 

Cl 

5.0 

1 


Cl 

3.7 

2 


C2 

-0.09 



Cl 

1.5 

3 


C2 

0.6 



Cl 

1.5 

4 


C2 

0.5 


DIFFUSION 

C S 

0.22 

5 


C S 

0.11 

6 

DISSIPATION 

C 

e 

0.15 



C £ 1 

1.44 

7 


Ce2 

1.92 
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Table 4.7 Comparison of calculated and measured locations 
of the stagnation points along the centerline in mm 


MODEL 

FORWARD 

REAR 

k-e 

31.4 

147.9 

ASM/MY 

50.8 

267 

ASM/R 

49.5 

242.7 

ASM/E 

49.7 

269 

RSM 

91 

275 

DATA 

38 

170 
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Table 4.8 


Various models used in flux transport 
and algebraic flux models 


MODEL 

MODEL NUMBER 

EQUATION 

NUMBER 

REDISTRIBUTION 

26 

(2.22) 


27 

(2.22) + (2.23) 

DIFFUSION 

28 

(2.20) 

VARIANCE 

29 

(2.25) 
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Table 4.9 Various flux transport and algebraic flux 


models 

(AFM) used in this 

investigation 

MODEL 

MODEL NUMBER 

EQUATION 

NUMBER 

AFM/MY 

30 

(2.36) + (model-26) 


31 

(2.36) + (model-27) 

AFM/ LAUNDER 

32 

(2.38) + (model-26) 


33 

(2.38) + (model-27) 

FLUX TRANSPORT 
MODEL 

34 

(2.8) + (model-26) 
+ (model-28) 


35 

(2.8) + (model-27) 
+ (model-28) 
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Table 4.10 


Values of constants used in scalar flux models 


MODEL 

CONSTANT 

ASSIGNED 

MODEL 



VALUES 

NUMBER 

REDISTRIBUTION 

c ie 

3.63 

26 


c ie 

C 2d 

3 . 0 
0.33 

27 

DIFFUSION 

C s* 

0.11 

28 

VARIANCE 

C Dd 2 

1.3 

29 



1.5 
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Figure 4.1 Sketch of inlet and test section A: Johnson and 
Bennett (1981), B: Roback and Johnson (1983) 
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Figure 4.2 Sketches of flow regions for swirling and 

non-swirling flow conditions A: Johnson and 
Bennett (1981) , B: Roback and Johnson (1983) 
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4.4 Comparison of predicted & measured axial 
velocity profile using k-* model 
(Johnson and Bennett (1981)) 
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Figure 4 . 6 Comparison of predicted & measured turbulent 
shear stress profile using k-< model 
(Johnson and Bennett (1981)) 
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Figure 4.7 Comparison of predicted & measured centerline 
axial velocity profile 
(Johnson and Bennett (1981)) 
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Figure 4 . 9 Comparison of measurements with k calculations 
using n.. from eq. (2.13) 

(Johnson and Bennett (1981)) 










Figure 4.11 Comparison of measurements with 

calculations using II.. from eq. (2.13) 
(Johnson and Bennett (1981)) 
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Figure 4 . 13 Comparison of measurements with y/w 2 

calculations using n.. from eq. ( 2 . 13 ) 
(Johnson and Bennett ( 1981 )) 
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Figure 4.14 Comparison of measurements with U and uv 
calculations using II.. * n n , + n.. 2 
(a) egs. (2 . 13 , 2 . 14) , J (b) egs. 13 , 2. 15) , 

(c) egs. (2 . 13 , 2 . 16) (Johnson and Bennett (1981)) 
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Figure 



.16 Comparison of the distributions of numerical 
viscosity with turbulent viscosities at 
different locations along combustor 
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Figure 4 . 17 Comparison of hybrid differencing scheme 

with quadratic upwind differencing scheme on 
mean axial velocity 
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Figure 4.18 Comparison of measurements with mean axial 
velocity calculations using two turbulent 
diffusion models (2.18,2.19) 


138 







Figure 4.21 Comparison of measurements with 

calculations using two turbulent diffusion 
models (2.18,2.19) 
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Figure 4.27 Comparison of measurements with k calculations 
(Johnson & Bennett (1981) , k-e, ASM, RSM) 
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Figure 4.31 Shear stress and normal stress production of 
turbulent kinetic energy 
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Figure 4.32 Secondary recirculation in the corner region 

using a high and a low Reynolds stress closure 








using a high and a low Reynolds 
number model (RSM) 
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Figure 4.35 Comparison of predicted & measured centerline 
axial velocity profile 
(Roback & Johnson (1983)) 
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Figure 4.36 Comparison of predicted & measured axial 
velocity profile using k-£ model 
(Roback & Johnson (1983)) 
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Figure 4.39 Comparison of measurements with mean 
tangential velocity calculations 
(Roback & Johnson (1983), ASM, k-e) 
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Figure 4.41 Comparison of measurements with 

calculations (Roback & Johnson (1983) , ASM) 
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Figure 4.46 Comparison of measurements with mean scalar 

calculations using AFM (Johnson & Bennett, 1981) 
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Figure 4.47 


Comparison of measurements with 
calculations using AFM (Johnson 


■J? 

& Bennett, 1981) 
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Figure 4.48 


Comparison of measurements with v6 
calculations using AFM (Johnson & Bennett, 1981) 
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Figure 4.50 Comparison of measurements with v$ 

calculations using second-moment closure 
(Johnson and Bennett , 1981) 
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Figure 4 . 52 Comparisons of measurements with v6 

calculations using AFM (Roback & Johnson, 1983) 
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AFM by Mel lor & Yamada, Pressure Redistribution Model by Monin 

AFM by Mellon & Yamada, Pressure Redistribution Model by Launder 

AFM by Launder & Gibson, Pressure Redistribution Model by Monin 

AFM by Launder & Gibson, Pressure Redistribution Model by Launder 

DATA (ROBACK AND JOHNSON, 1983) 


Figure 4.53 Comparisons of measurements with v0 

calculations using AFM (Roback & Johnson, 1983) 
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CHAPTER 5 

REACTING FLOW CALCULATIONS 

This chapter presents the results of the reacting flow 
calculations using the constant-density turbulence models 
recommended in the previous chapter. There is still the 
basic question of whether empirical equations or models 
developed for constant density, non-reacting flows can be 
applied without modification to turbulent combustion. Models 
are generally required for quantities such as the turbulent 
eddy viscosity and the scalar dissipation function which 
plays an important part in turbulent chemistry (Bray, 1980) . 
There is a widely held belief that the use of Favre 
averaging automatically takes account of all the effects of 
density fluctuations and chemical reactions, and, therefore, 
allows empiricism from constant-density flows to be 
exploited without change. While this may indeed be true in 
some cases, it can be justified only by actual comparison 
with experiments. 

Another fundamental question in the modelling of 
reacting flows is concerned with the effects of turbulence- 
induced fluctuations on time-averaged reaction rate. An 
essentail feature of many of the practical problems is that 
the rate of chemical reaction is limited by the rate of 
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mixing of reactants. For non-premixed system, and in the 
limit of fast chemistry, the molecular species are 
instantaneously related to the conserved scalar and the 
statistics of all thermodynamics variables can be determined 
from a sufficient knowledge of the statistics of that scalar 
(Bilger, 1976) . As a result, the need to evaluate mean 
reaction rate is obviated. However, the estimation of 
unburnt fuel and the formation of CO require consideration 
of finite-rate reactions. On the other hand, premixed 
combustion, in contrast to non-premixed combustion, requires 
the evaluation of mean reaction rate because the mixing is 
accomplished before reaction begins. Based on these 
differences, the results to be presented are divided into 
non-premixed and premixed combustion. Under non-premixed 
combustion, both results of fast and finite rate chemistry 
models are discussed for flows with and without swirl. 

In view of turbulence/chemistry modelling, it would be 
very usefuel if isothermal and reacting flow experiments 
could be carried out in the same apparatus so that a natural 
progression in understanding could be followed. Generally, 
because of experimental problems, such progression has not 
been done. Therefore, effort is made to select the test 
cases which are similar to the previous isothermal flow 
cases in terms of the geometry, working fluids, boundary 
conditions and measurements technique. In addition, all data 
sets have to meet the criteria for acceptable measurements 
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in terms of adequate experimental facility, appropriate 
instrumentation and agreement with generally accepted flow 
trends. In view of the above considerations the following 
three experiments are chosen as the test cases: 

1. Lewis and Smoot (1981) 

The configuration of this experiment is very close to the 
configuration of Johnson and Bennet (1981) . 

2. Brum and Samuelsen (1982) 

This experiment is used to simulate the dilute swirl- 
stabilized combustor. Since molecular weights of C0 2 and 
C 3 H 8 are identical, therefore, by comparing with the non- 
reacting results, the effect of heat release on flow field 
can be assessed. 

3. McDannel et al. (1982) 

This experiment is used to simulate an axisymmetric opposed 
reacting jet combustor. 

Again, the shortcoming of these experiments is the lack of 
well-defined inlet boundary conditions. Nevertheles, the 
calculations are quite beneficial because they help to 
illustrate the differences between the various models and 
their ability to calculate reacting flows. 
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5.1 NON-PREMIXED COMBUSTION 


5.1.1 Coaxial Non-swirling Jets 

In the this study, the geometry of Lewis and Smoot 
(1981), simulating an industrial furnance, is selected. In 
this experiment, coaxial streams of fuel (town gas ) and air 
are injected into a suddenly-expanded combustion chamber 
(Figure 5.1). The flame is stabilized at the dividing lip 
between the two streams. Measurements have been made of the 
time mean mixture fraction and species concentration. The 
parameters and test conditions for this combustor are 
summarized in Table 5.1. 

The computations for this case are made assumption a 
two-dimensional formulation and a standard k-e model. This 
model is previously recommended for non-swirling, constant- 
density flows. Two modelling approaches have been used for 
the combustion process. At one condition, chemical kinetics 
are assumed to be rate controlling and at another, turbulent 
mixing is treated in detail, but infinite rate chemistry is 
assumed. A non-uniform grid of 61x47 is used in the 
computations. Uniform axial velocity profiles are prescribed 
for the fuel and air stream, respectively. The inlet 
turbulence intensities for air and fuel are given in the 
measurements as 6% and the length scales are assumed to be 
5.7 mm for the air and 1.6 mm for the fuel jet, 
respectively. Along the adiabatic walls, standard wall 
function treatment is employed. The fuel mixture fraction is 
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set equal to one in the fuel stream and zero elsewhere. 
Because of the uniform distribution of mixture fraction in 
each of the streams, the fluctuation values are zero in the 
inlet plane. 

The predictions of the mixture fraction, unburnt fuel 
and carbon monoxide using fast and finite-rate chemistry 
models are presented in Figure 5.2. The rate constants used 
are given in Table 5.2. It can be seen that the results are 
qualitatively correct. Both models over-predict the mixture 
fraction near the centerline in the developing region; 
however, the finite-rate chemistry model seems to be able to 
reproduce the physics better than the fast-chemistry model. 
A turbulent Prandtl number of 0.9 is used to obtain the 
results shown. In the fully developed region, the agreement 
is quite good. This suggests that a good prediction of 
mixture fraction in the developing region requires a lower 
turbulent Prandtl number. In the region near the inlet, the 
convection is dominant and the change of the Prandtl number 
does not have any effect on the prediction of mean mixture 
fraction. In the developed region, since the absolute levels 
of mixture fraction are low, the difference caused by 
changing the Prandtl number will not be significant. The 
middle zone is the region of very steep variation in mixture 
fraction and a lower Prandtl number will undoubtedly improve 
the prediction. Comparison of predicted and measured values 
of unburnt fuel is similar to that of mixture fraction. The 
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prediction is qualitatively correct at all planes; however, 
the model has over-estimated the unburnt fuel in the middle 
region indicating the slow reaction rate. The predicted CO 
levels are fair in most of the locations but beyond x =785 
mm, they are significantly lower than the data. The 
predicted results shown in Figure 5.2 are based on the 
assumption of equal diffusivity for all species, enthalpy, 
and the mixutre fraction. It is obvious that the values of 
the turbulent Prandtl number are not necessarily the same 
for all species and may not even be uniform over the whole 
flow field. Therefore, further investigation regarding the 
scalar transport model is required. 

The contour plots of temperature, mixture fraction, 
unburnt fuel, carbon dioxide and carbon monoxide inside the 
combustor are shown in Figures 5. 3-5. 7. The flame zone is 
clearly delineated by the temperature contours. The 
comparison clearly shows the big difference between the two 
combustion models, and illustrates the strong interactions 
between flow and chemical reactions for the finite-rate 
chemistry model and the relatively weak interactions shown 
in the fast-chemistry model (Figure 5.3). Furthermore, a 
sharp temperature gradient seen near the edge of the flame 
is a direct result of an insufficient diffusion along the 
radial direction. No doubt, a smaller Prandtl number will 
help. However, a more realistic approach would be to abandon 
the constant Prandtl number assumption and proceed to 
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evaluate the turbulent fluxes by their own transport 
equations. The carbon monoxide profile along the centerline 
behaves correctly, attaining a maximum when reaction (2.63) 
and (2.64) are equally competing, and subsequently falling 
off as the combustion products tend to their equilibrium. In 
general, the finite-rate chemistry model seems to give more 
realistic results and could be easily adopted for non- 
swirling combustor flow calculations. 

The cold flow results are also shown in Figures 5. 8-5.9 
for comparison, and serve to show the effect of heat release 
on the mean and turbulence field. The recirculation zone in 
the case of reaction is: (1) more intense (higher negative 
velocities) , and (2) more compact (shorter) . Expansion 
effect due to heat release is amply demonstrated by the 
substantially lower mean axial velocity (U) in the developed 
region of the flow. On the other hand, turbulence activities 
as measured by k and uv are greatly reduced by heat release. 
One surprising result of the present study is given by the 
calculation of the fast chemistry model. Essentially, a 
similar behavior of U, k and uv are obtained in the reacting 
zone in comparison with the isothermal flow calculations. 
One reason for this could be due in part to unrealistic 
simulation of heat release effects on density fluctuations 
in the transport equations. It also leads to prediction of 
double peaks behavior for the mean axial velocity 
immediately after the inlet plane which is not physically 
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correct. Such a behavior might be attributed to the 
predicted scalar distribution near the inlet plane, 
especially in the wake region created at the dividing lip 
between the two streams. In a mixing controlled combustiom 
model, the hydrodynamic field is highly sensitive to the 
mixture fraction prediction. It is the only parameter that 


determines 

the 

amount 

of heat 

release . 

An 

incorrect 

prediction 

of 

mixture 

fraction 

results 

in 

incorrect 


temperature and density fields. Since density is strongly 
coupled to the flow field, an incorrect density changes the 
flow field completely. On the other hand, in finite rate 
chemistry model, the chemical heat release is determined by 
transport of total enthalpy and unburnt fuel which is a more 
realistic approach. It must be stated that neither models 
can provide an entirely satisfactory description of the 
reaction zone, however, the two-step reaction scheme has 
more flexibility and greater potential for application in 
gas turbine combustors. It needs to be further validated 
with simple flames to establish rate constants. In addition 
to Favre-averaging technique, the turbulence model has to 
include the heat release effect on the turbulence field in 
the transport equations of turbulent kinetic energy and 
scalar fluctuation too (Dibble et al. 1985 ) . 



5.1.2 The Dilute Swirl-Stabilized Model Combustor 

The flow geometry in this case is identical to that 
shown in Figure 4.3. The propane-fueled combustor is 
operated at atmospheric pressure with an overall equivalence 
ratio of 0.2 and a bulk reference velocity of 15 m/s. The 
fuel velocity based on the nozzle area and mass flow is 28.1 
m/s which is high and has a tremendous effect on the 
prediction of the central recirculation zone. 

Computations are made over a 61x47 non-uniform grid 
using a standard k-e and an equilibrium ASM model. The rate 
constants used in this calculation are given in Table 5.2 
and Prandtl/Schmidt number is taken to be 0.9. Predicted 
velocity results from the k-6 and the equilibrium ASM models 
for isothermal and reacting flows are shown in Figures 5.10- 
5.12. The effect of swirl is pronounced. It produces a large 
recirculating region near the axis that sweeps back burnt 
gas products to ignite the incoming reactants. A transition 
occurred in the form of the recirculation zone which becomes 
shorter and wider. Beyond the recirculation region, the mean 
axial velocity profiles predicted by the ASM/E are in better 
agreement with data than those predicted by the k-e model. 
However, both models failed to predict the correct location 
and strength of the recirculation region along the 
centerline. This result is attributed partly to the high 
vertical component of inlet fuel velocity used which tends 
to push the central recirculation region towards the 
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combustor wall. The predicted maximum flow reversal velocity 
is approximately -5.75 m/s, while the corresponding measured 
value is about -7.95 m/s. Furthermore, the location of the 
maximum negative velocity is predicted closer to the swirler 
exit compared to the data. In the developing region, the 
ASM/E has under-predicted the turbulent intensity, 
especially in the central recirculation zone; however, in 
the fully developed region the results are relatively better 
(Figure 5.11). It is also observed that, at downstream 
locations along the centerline, the experimental value of 
is greater than . These excessive levels as mentioned 
by Brum and Samuelsen (1982) are probably produced by a 
spiralling action of the vortex center about the stationary 
laser velocimeter probe which can be considered as an 
experimental error. The discrepancy between prediction and 
measurement can be attributed to: (1) inaccurate inlet 
conditions, (2) the turbulence model constants, (3) 
experimental error and (4) the effect of heat release on the 
turbulence field which needs to be accounted for ( Dibble et 
al. , 1985). 

The effect of reaction is observed to double the mean 
axial velocity and raise the turbulence intensity in the 
reacting zone (Figure 5.10). However, the increase of 
tangential velocity near the centerline, which also stems 
from heat release, has not been correctly predicted. Unlike 
the previous test case, the fast and finite-rate chemistry 
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results are essentially identical. The reason is the 
enhanced mixing due to swirl. This increased mixing activity 
promotes chemical reactions between fuel and oxidant, and 
the net result is one of very fast reaction rates. The 
recirculation zone in the case of reaction is wider (16 mm) 
and the length of the recirculation region is shorter (38 
mm) . These changes are the results of gases expanding in the 
shear layer at the boundary of the recirculation zone. This 
creates a region of low pressure which tends to draw 
swirling air towards the centerline and close the 
recirculation zone (Brum and Samuelsen , 1982) . 

Finally, contours of unburnt fuel, temperature, mass 
fraction and carbon dioxide are shown in Figures 5.13-5.16. 
There is a strong gradient in fuel concentration in the 
radial direction close to the inlet, while the profile 
becomes more uniform as the outlet is approached. The 
temperature and mixture fraction distributions found by fast 
and finite rate chemistry models are very similar. However, 
the fast chemistry model shows a longer reacting zone. Both 
models predict the sharp temperature gradients near the 
perimeter where dilute air enters. This apparent suppression 
of radial transport could be due in part to the use of an 
isotropic diffusion model, unrealistic Prandtl/Schmidt 
number, and unrealistic simulation of heat release effects 
on density fluctuations in the transport equations. As 
mentioned before, the Prandtl number is not a significant 
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factor in the fully developed region because the flow is 
convection dominated and the axial velocity is quite 
uniform. However, in the developing region, the isotropic 
diffusion model is not valid and more realistic approach is 
to evaluate the turbulent fluxes by their respective 
transport equations. 

5.2 PREMIXED COMBUSTION 

Another test case selected for evaluating the kinetic 
scheme is the recirculating flow in an axisymmetric, 
opposed-reacting- jet combustor shown schematically in Figure 
(5.17). The incoming mainstream of premixed propane and air 
(U m =7.5 m/s) is opposed by a high-velocity (Uj=135 m/s) 
premixed jet positioned along the longitudinal axis with an 
equivalence ratio of one. The jet creates a zone of 
recirculating flow necessary to stabilize the reaction. 
Computations are made over a 49x25 non-uniform grid using a 
k-£ model and the results compared to species concentration 
and temperature measurements reported by McDannel et al. 
(1982) and velocity data provided by Samuelsen (1986) . The 
inlet turbulence intensity is assumed to be uniform with a 
value of 0.005 and the inlet length scale is assumed to be 
constant with a value of 2.5 mm. Adiabatic boundary 
conditions along the walls are also applied. Again, the rate 
constants listed in Table 5.2 are used. 
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The difficulty of modelling jets in a co-flowing stream 
is amply demonstrated by the work of Launder et al. (1973). 
Their results show that, of the four models, ranging from 
mixing-length to two-equation models, used to calculate jets 
in a co-flowing stream, none is capable of predicting the 
decay of the jet centerline velocity correctly. So and Hwang 
(1986) point out that there are certain inherent 
deficiencies in the models and these lead to a slower growth 
for the jet. A model for the developed region of the jet is 
put forward by So and Hwang (1986) and their calculations 
in good agreement with measurements. For opposed— jets in 
a moving stream, the modelling difficulties are further 
compounded by the presence of a stagnation point, and the 
very rapid decay of the jet centerline velocity. In other 
words, the models may not be able to mimic the highly 
dissipative phenomenon occurring in the region downstream of 
the jet exit. In view of this, the calculated isothermal 
flow field of the opposed-jet experiment could not be 
expected to be correct. Therefore, the present comparison 
with measurements should be judged bearing in mind the 
inadequacy of the turbulence model and the inaccuracy of the 
isothermal flow field. 

A comparison of the calculated jet centerline velocity 
decay with the hot- and cold-flow measurements is shown in 
Figure 5.18. As expected, the models under-predict the 
velocity decay for cold flow and the discrepancy becomes 
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more pronounce for the hot flow. The measurements show that 
the jet decays immediately downstream of the jet exit, while 
the models give a very small potential core region for the 
cold-flow jet and a much larger potential core for the hot- 
flow jet. The location of the stagnation point is also 
incorrectly predicted, and the result shows that the 
predicted hot-flow stagnation point is in closer agreement 
with the cold-flow measurement. Experimental measurements 
show that the hot-flow centerline velocity is always higher 
than the cold-flow centerline velocity. This is a 
consequence of flow expansion due to heat release from 
chemical reaction. The models cannot reproduce this behavior 
correctly because they over-predict the centerline velocity 
decay in the region 1.7<x/D<2.6 for the hot-flow and under- 
predict the decay in the same region for the cold-flow. As 
for the region between the stagnation point and the jet 
exit, the k-e model completely fails to reproduce the flow 
behavior. In the region downstream of the stagnation point, 
the agreement is much improved. Therefore, this points to 
another deficiency of the k — e model; namely, it cannot 
properly account for the rapid dissipative phenomenon seen 
in the present test case. 

The predicted mean axial velocity profiles for hot and 
cold flows are presented in Figure 5.19. Shown for 
comparison are the hot— and cold-flow measurements obtained 
by Samuelsen (1986) at one x location. The expansion effect 
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due to heat release is clearly apparent. The calculated 
recirculation zone in the case of reaction is (1) radially 
wider, (2) axially shorter and (3) more intense, which means 
higher negative velocities. This is in stark contrast to the 
measured profiles which show a radially wider recirculation 
zone for the cold flow. Reason for this could be traced to 
the inadequacy of the k-t model and perhaps, partially, to 
the combustion model. The effect of reaction is also to 
increase the turbulent kinetic energy and shear stress 
considerably (Figures 5.19-5.20). In both cases, the 
stagnation point is incorrectly located due to deficiencies 
in the underlying turbulence model. 

Not only does the turbulence model affect the flow 
field calculation, it also affects the temperature and 
species concentration field. A comparison of the temperature 
and carbon monoxide (CO) distribution for stoichiometric 
(<£=1) combustion of propane and air is shown in Figures 
5.21-5.22. The flame region is indicated clearly by the 
temperature contours, but local agreement is not attained. 
As discussed before, the stagnation point error is clearly 
evident and so are the under-predicted axial turbulent 
exchange across the stagnation point and the deficient 
radial exchange. 
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5.3 CONCLUSIONS 

The main conclusions emerged from this study can be 
summarized as follows. 

1. Two-step reaction scheme performed better than fast 
chemistry model in predicting mean mixture fraction. 

2. Satisfactory prediction of mean mixture fraction can be 
obtained using a turbulent Prandtl number that is constant 
over the flow field. 

3 . Finite-rate chemistry models are preferred over fast- 
chemistry models for determining the combustion effects in 
combustors. However, it needs to be further tested with 
simple cases to establish rate constants so that major 
species can be accurately predicted. 

4. The constant-density k-e turbulence model provides a 
satisfactory representation of the aerodynamics in most 
practical combustor flows, except in the case of jet- 
stabilized combustor flow. The reason is the inability of 
the k-e model to replicate the highly dissipative phenomenon 
found in such flows. 

5. Some improvement in complex swirling flow predictions 
could be obtained by using an algebraic stress model. 

6. The combustion model must provide a realistic 
representation of heat release behavior which has the 
substantive effect on the flow structure. 
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7. In addition to Favre-averaging, the turbulence model 
should include some heat release effect on the turbulence 
field. This could be accounted for in the transport 
equations of turbulent kinetic energy and its dissipation 
and the scalar fluctuations. 
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TABLE 5.1 Combustor parameters and test conditions 


Primary Gas 


Temperature (K) 

Mass flow rate (Kg/s) 
Velocity (m/s) 
Composition (molar %) : 


286 

0.0031 

21.7 


ch 4 

88 . 53 

C 2 H 6 

7.44% 


2.55% 

CO-, 

1.39% 

H 2 

0.09% 


Secondary Air 


Temperature (K) 

Mass flow rate (Kg/s) 
Velocity (m/s) 
Composition (molar %) : 


589 

0.0362 
34 . 7 


2 

2 

Ar 


78 . 3 
20.8 
0.9 
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TABLE 5.2 Model constants for mean reaction rates 


Town Gas Propane 


C R1 

" 

3 . 

3 . 


C R2 

= 

4 . 

4. 


H 

* 

= 

3.3 X 10 14 

5.57 X 

10 

K 2 

= 

00 

o 

H 

X 

o 

• 

kd 

5.42 x 

10 

El/R 

= 

27000 

15104 


E2/R 

= 

12500 

15098 


a 

= 

1.5 

1.75 


b 

= 

0.5 

0.1 


c 

= 

1.0 

1.65 


a' 

= 

2 . 0 

2.0 


b' 

=r 

1.0 

1.0 


c ' 

= 

1.0 

0.5 


d' 

— 

0.0 

0.5 
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Sketch of the combustor (Lewis & Smoot, 1981) 
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Figure 5.2 A Comparison of the calculated and measured 
mixture fraction, fuel and CO 
(Lewis & Smoot ,1981) 
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Figure 5.3 Contour plots of temperature distribution 
inside the combustor (Lewis & Smoot, 1981) 
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Figure 5.4 Contour pi 
(Lewis & S 
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Figure 5.5 Contour plots of unbumt fuel 
(Lewis £ Smoot, 1981) 
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Figure 5 . 6 Contour plots of carbon dioxide 
(Lewis & Smoot, 1981) 
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Figure 5.7 Contour plots of carbon monoxide 
(Lewis & Smoot, 1981) 


204 


U (m/s) 



k (m/s)^ 


Figure 5.8 Effect of heat release on the calculated mean 
and turbulent flow fields 
(Lewis & Smoot, 1981) 
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Figure 5.12 Comparison of measurements with uw calculations 
(Brum & Samuel sen (1982) , ASM) 
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Figure 5.13 


Contour plots of unburnt fuel 
(Brum & Samuelsen ,1982) 
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Figure 5.14 Contour plots of temperature distribution 

inside the combustor (Brum & Samuelsen ,1982) 
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Figure 5.15 Contour plots of mixture fraction 
(Brum & Samuelsen ,1982) 
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Figure 5.16 Contour plots of C0 2 (Brum & Samuelsen ,1982) 
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Figure 5.18 Opposed jet centerline velocity decay for hot 
and cold flow and their comparison with 
calculations (McDannel et al. ,1982) 
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Figure 5.19 A comparison of the axial velocity 

profiles at different axial locations 
(McDannel et al. , 1982) 
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Figure 5.21 Comparison of the calculated and measured 
temperature contours inside combustor. 
(McDannel et al. , 1982) 
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Figure 5.22 Comparison of the calculated and measured CO 
contours inside combustor. 

(McDannel et al. , 1982) 
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CHAPTER 6 


CONCLUDING REMARKS AND RECOMMENDATIONS 


Important specific conclusions for each of the sections 
in chapters 4 and 5 are provided at the end of each section. 
This chapter, therefore, presents general conclusions of 
this study and makes recommendations for further work. 


6.1 CONCLUSIONS 

1. As far as the mean field prediction is concerned, k-e 
closure performs just as well as any ASM or RSM in the 
calculation of non-swirling combustor flows. 

2 . The k-e model gives good correlation for the developing 
region of complex swirling flows. However, for the far-field 
region the ASM provides a better prediction. 

3. For swirling combustor flows, the ASM's and the full 
Reynolds stress models do a better job of predicting 
tangential velocities, while the k-e model gives a good 
description of the centerline recirculation zone although 
the predictions are subject to uncertainties from the inlet 
boundary conditions. 
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4. Low-Reynolds-number model provides a better estimate of 
the mean and turbulence quantities in the near wall region. 
This model also predicts the corner recirculation zone which 
is observed experimentally and missed completely by all 
high-Reynolds-number models. 

5. Models employing transport equations for the individual 
turbulent stress and fluxes components simulate the 
turbulent processes more realistically and are therefore 
potentially more general compared to the simpler models. 
However, they are not thoroughly tested and are 
computationally more expensive. Hence, at the present state 
of development, they are not very suitable for practical 
applications. They are important, however, as a starting 
point for deriving algebraic expressions for the turbulent 
stresses and fluxes. It seems that such expressions used in 
conjunction with the k and e equations are sufficient for 
most engineering problems. 

6. It should be noted that the effectiveness of turbulence 
model predictions could be obscured to some extent by 

competing factors ' boundary conditions, oscillatory 

phenomena and numerical diffusion. A significant 
contribution from any of the aforementioned factors tends to 
invalidate conclusions regarding the superiority or 
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inferiority of a given turbulence model. Numerical diffusion 
is a complex function of mesh size and cell aspect ratio. 
Inlet and boundary conditions are also of importance in 
strongly swirling flows. The turbulence model cannot 
compensate for inadequacy in this area. It is apparent that 
turbulence model validation must be preceeded by (1) 
appropriately configured and detailed experimental case 
studies, and (2) elimination of false diffusion 
considerations . 

7. Favre-averaging technique is a reasonable approach for 
isothermal, variable density flows. However, for reacting 
flows, the turbulence model should also include the effects 
of chemical heat release on the Reynolds stress/flux 
components . 

8. Two-step reaction scheme shows promise for application in 
gas turbine combustors and is preferred over fast-chemistry 
model for determining the combustion effects in combustors. 
However, they have to be further validated with simple 
flames to establish model constants and rate constants, so 
that the major species can be accurately predicted. 
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6.2 RECOMMENDATIONS 


1. The low-Reynolds-number closure is found to provide a 
fritter estimate of the mean and the Reynolds stresses near 
the wall region. It is appropriate to apply a similar 
approach to develop a low-Reynolds-number model for the 
scalar transport equations too. 

2. Although some of the models described in this review and 
in particular the k-e model, have shown to work well in many 
situations, there is much room for further developments. The 
e-equation in its present form appears not to be 

iciently universal and should be improved. As observed 
by many, this equation is the Achilles heel for most models. 
Ideas to use several length-scale equations for different 
directions or different processes are promising (Hanjalic et 
al. 1979) and should be developed further. 


3 . The model 

assumptions for the 

pressure-strain 

and 

pressure-scalr 

gradient 

correlations 

are also 

not 

very 

satisfactory 

and need 

improvement. 

Proposals 

for 

the 


behavior of , 1 in inhomogeneous flows have only gone 
further than that of Rotta (1951) by including further terms 
in a series expansion about the isotropic homogeneous state 
(Lumely and Khaheh-Nouri , 1974). However, optimization of 
the coefficients of the terms in the expansion on the basis 
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of available experimental data is a very difficult task 

indeed. It seems unlikely that any serious proposal for 

n.. will be made in the near future and that emphasis 

3-D 1 1 

should be placed on developing a better approximation for 
n ij,2‘ 

4. The derivation and validation of higher-order closure 
schemes holds the greatest potential for turbulence model 
improvement for strongly swirling flows. 

5. The difficulty in clearly demonstrating calculations free 
of numerical errors has restricted the testing of higher- 
order closures in recirculating flows. Efforts to find a 
stable and higher order (order of terms retained in an 
equivalent Taylor series expansion) differencing scheme that 
can eleminate numerical diffusion should continue. This is 
especilly important in the case of reacting flows because of 
the coupled non-linearities which exist between the chemical 
and fluid mechanical processes. 

6. The correct solution of the potential core region is very 
important in analyzing the chemical heat release effects at 
the early stages of the reaction and obtaining the correct 
behavior of axial variations of turbulent velocity/scalar 
fluctuations in the diffusion flames. 
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7. An intensive submodel validation and development efforts, 
especially for the algebraic stress/flux model, two-step and 
four-step kinetic schemes, and two-step kinetic scheme in 
conjunction with PDF approach should be continued. 

8. Unfortunately, in many instances there is a lack of 
quality data relevant to gas turbine combustion. Many 
modelling assumptions are similar to the constant— density , 
Reynods-stress closure. Therefore, further experiments with 
more emphasis on turbulent scalar fluxes and density 
correlations are needed to support or to improve these 
assumptions. In addition, more experiments with different 
fuels are required to assess the idealized density-mixture 
fraction relation and its application to turbulence 
modelling of diffusion flames. 
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APPENDIX A 

TURBULENT FLOW EQUATIONS FOR THE k-e MODEL 

The transport equations presented in the previous 
sections reduce in axisymmetric coordinates (x,r) to the 
following 

CONTINUITY EQUATION : 
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APPENDIX B 


TURBULENT FLOW EQUATIONS FOR THE ASM MODEL 


The governing equations under the gradient diffusion 
assumption in the cylindrical coordinate system are reduced 
to the following 

CONTINUITY EQUATION : 
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APPENDIX C 


ALGEBRAIC-STRESS CLOSURE IN AXISYMMETRIC COORDINATES (x,r) 

The Reynolds stresses are obtained by the following 
algebraic expression 


P. % bS. . t 

1 J 3 l j 


au„ 


C. r(u.u.- = (. .k) — a ( P . =6 . .P, - xu.u. ) 

1 k v i j 3 ij ' ij 3 ij k 3 l j dx ffl ' 

(Cl) 

= 0 . 


au_ 


3‘i/k - 3 u i u j 57^> - rkS ij- k u i u j< P k~° 


Constants A and B are equal to 1 for Rodi’s model. Mellor 
and Yamada’s model is obtained by assigining A=0 and B=P k /c 
and for the equilibrium model A=0 and B=1 . In the 
cylindrical coordinates system, production of the Reynolds 
stresses and the mean strain rate are expressed as 
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APPENDIX D 


ALGEBRAIC FLUX MODEL IN AXISYMMETRIC COORDIANRES (x,r) 


The turbulent fluxes are obtained by the following 
algebraic expressions 


A^ (p k-O= 0 
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_ * — ae ae . 1W — - av j — - av w. 

If k 3x dr v 2f ' x dx Sr r 


+ A ^< P k-‘>=° 


( D2 ) 


. t — — ae — ae . , . aw . — - aw — v w. 

If k dx dr v 2f dx dr r r' 


+< 5sr< p k-‘>=° 


(03) 


Constant A is equal to 0 for Mellor and Yammada’s model; and 
is equal to 1 for Launder and Gibson ’ model . 
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APPENDIX F 

THE SCALAR FLUX TRANSPORT CLOSURE 



0u0 , —5- du 6 v 

uv + v 2 ) 


0x 


6r 


i-ifei 


k — dv 6 — v 2 

rpC . - ( 2 uv - w# - 

r sB t dx r 


— w 0 \ 1 . w —r - - r -r av a - s- av n — w — a© -.ae' 

vw — ) J + r ” e p - p L u * 33 T + ve 3 F _ 2w<? ? +uv 57 +v 3 F 


- c l# k v * 


Equation for wtf : 


1 d ( rpUwg ) 6 ( rpVwfl ) 

r dx 0r 


(F 3 ) 


l-Mfc 


(rpC 


k —5- 0 w 0 
- u 2 

s 6 t ox 


) + 


0 . - k — 0W0 . 

aF (r ' ,c s * I v eF - J 


J-fe[ 


- k , — 0w0 v 0 

' c .» 7 (uv SF~ + u ” — + 


du 6 0u0 N 

uw + vw — ) 


0X 


0r 




- k , 0w0 v B 

r ' c .* 7 (uv 57- - v “ — + 


UW 


dve ^ — dve w* .] - f — a© ^ — a© ^ — aw 

dx dr r J dx dr dx 


vtf jp- + w* - + v$ - 1 - C. p £ wff 
dr r r J 10 ^ k 


APPENDIX G 


BOUNDARY CONDITIONS FOR HIGH REYNOLDS -NUMBER MODELS 

Differential transport equations for various turbulence 
quantities were introduced in chapter two. There exist an 
infinite number fields that satisfy the same set of 
equations, and only the boundary conditions distinguish one 
type of flow from another. Thus consideration of boundary 
conditions is very important. Most turbulence models are 
devised for high Reynolds number flows. However, in the 
vicinity of the solid boundaries the low Reynolds number 
effects become significant and must be accounted for. This 
can be performed either by solving the low Reynolds number 
version of the transport equations or by developing wall 
functions that introduce these effects into the existing 
high Reynolds number models. The first option predicts the 
near wall region better but it needs vast amount of computer 
time. In the high Reynolds models, it is assumed that the 
distribution of mean velocity and turbulence quantities 
throughout the main part of the flow is weakly dependent on 
the model of the turbulent transport in the immediate wall 
vicinity and the conventional logarithmic law which is based 

| 

on the local equilibrium is applied. Details of this sort of 
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bounadary condition for high Reynolds number models are 
given in the following sub-sections. 


G.l NEAR-WALL VELOCITY PROFILE 

In general, near the wall, one dimensional couette flow 
analysis is made. In this region shear stress is assumed to 
be constant. This condition is true only for an impermeable 
wall, with negligible streamwise pressure gradient. The wall 
region is made up of three zones: The viscous sublayer, 
transition or buffer zone, and inertial sublayer. Our 
approach is to dispose of the buffer layer by defining a 
point y + =11.63 below which the flow is assumed to be purely 
viscous and above which it is purely turbulent. The law of 
the wall may be expressed as 


y 11.63 


y > 11.63 


where 


+ 

U = 


+ 

Y 


+ 1 * 4 - 
U = ± £n(E y) 


(G 1) 
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P 


k , = Van Karman conastant = .4187 
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E = Roughness Parameter = 9.973 


G. 2 NEAR-WALL TURBULENCE KINETIC ENERGYY BUDGET AND 
DISSIPATION RATE 


The near-wall kinetic energy levels are obtained from 
the solution of the k transport equation. The convection and 
diffusion terms are treated by hybrid differencing while the 
production and dissipation terms need to be changed to 
include the near wall effects. The approach is based on one 
dimensional, constant shear stress couette layer. In the 
inertial sublayer and in the absence of buoyancy effects, 
the local rate of production of turbulence is balanced by 
the viscous dissipation rate e which, together with eddy 
viscosity and the fact that the shear stress is 
approximately equall to the wall shear stress r w , leads to 

2 



(G2) 


From equation (G2) the surface shear stress may be found as 


.25 

, p k \/k U C 

r s *— F p 

£n(E y) 


(G3) 
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The mean production rate of k per unit volume can then be 
found as 

MEAN PRODUCTION RATE = r s U p /y p (G4) 

where y p being the distance of the near wall grid point from 
the the wall and U p is the velocity at that point. 

In wall-flows, unlike turbulent kinetic energy which 
falls to zero at the wall , e reaches its highest value at 
the wall. With local equilibrium assumption 


and dU/dy from (G3) , there results the following boundary 
condition for e: 



(G5) 


G . 3 NEAR-WALL SCALAR PROFILE 


For an impermeable wall with zero streamwise pressure 
gradient, the scalar flux at pont P near the wall can be 
written as 



6 P ViZsP) 


(G6) 
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where 


e + = o 9 ' t (u + + PF) (G7) 

The quantity PF in Eq. (G7) is a following Van Driest 
Function: 

— _ 5 —PS 

PF = ( . 257r/Sin(jr/ 4 ) ) (/c/26) (-1 + a/<r t ) (a/a t ) * (G8) 

where a and cr t are the the molecular and turbulent prndtl 
number for scalar 

G . 4 SYMMETRY CONDITIONS 

At symmetry plane, the normal gradients are zero for 
all quantities with symmetrical behavior such as scalar 
quantities and velocity components parallel to the symmetry 
plane or line 

G . 5 INLET CONDITIONS 

For flows without swirl, the radial component of 
velocity is set to zero and axial velocity can be specified 
either by a fully developed pipe flow condition or by 
information from expeimental data. In swirl flows, howevere, 
an accurate specification of velocities is required for a 
better comaprison with measurement. The turbulent kinetic 
energy profile either is calculated by the mixing length 
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hypothesis or is determined from the experimental data. The 
inlet dissipation rate is calculated from the following 
expression 



where i is the turbulent length scale. 
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